Deterministic Algorithms for Four-Dimensional Imaging in Colocated MIMO OFDM-Based Radar Systems

In this manuscript, the problem of detecting multiple targets and jointly estimating their spatial coordinates (namely, the range, the Doppler and the direction of arrival of their electromagnetic echoes) in a colocated multiple-input multiple-output radar system employing orthogonal frequency division multiplexing is investigated. It is well known its optimal solution, namely the joint maximum likelihood estimator of an unknown number of targets, is unfeasible because of its huge computational complexity. Moreover, until now, sub-optimal solutions have not been proposed in the technical literature. In this manuscript a novel approach to the development of reduced complexity solutions is illustrated. It is based on the idea of separating angle estimation from range-Doppler estimation, and of exploiting known algorithms for solving these two sub-problems. A detailed analysis of the accuracy and complexity of various detection and estimation methods based on this approach is provided. Our numerical results evidence that one of these methods is able to approach optimal performance in the maximum likelihood sense with a limited computational effort in different scenarios.

INDEX TERMS Dual-function radar-communication, four-dimensional radar imaging, frequency estimation, multiple-input multiple-output radar, orthogonal frequency division multiplexing, radar processing.

I. INTRODUCTION
W IRELESS communication and radar sensing have been advancing independently for many years, even though they share various similarities in terms of both signal processing and system architecture. This consideration and the problem of radio spectrum scarcity have motivated the investigation of a new class of wireless systems, able to accomplish sensing and communication jointly. Various recent research activities in this field have evidenced that joint communication and sensing (JCAS) systems can provide significant advantages in terms of device size, power consumption, cost, and spectral efficiency compared to traditional systems in various applications [1]. Different approaches to their development are currently being investigated [2], [3], [4], [5], [6]. In this manuscript we adopt a communication-centric approach; this means that radar sensing represents an add-on to the considered wireless communication system. Moreover, we assume that orthogonal frequency division multiplexing (OFDM) is employed for both communication and sensing and that our system is equipped with both transmit (TX) and receive (RX) arrays (i.e., it is of multiple-input multiple-output, MIMO, type).
The OFDM modulation format has been adopted in various wireless communication standards, due to its robustness to multipath fading and to its relatively simple synchronization [5]; in addition, its use in MIMO communication systems has been widely investigated (e.g., see [7], [8] and references therein). A wide technical literature on the use of OFDM for radar sensing refers mainly to single input, single output (SISO) systems [2], [9], [10], [11], [12], [13], [14], [15], [16]. For this class of systems, various direct and indirect sensing methods for target detection and estimation have been proposed. Generally speaking, direct sensing methods extract target information from the received signal without compensating for the effect of the data payload it conveys (e.g., see [1], [11], [17]) and typically exploit computationally intensive compressed sensing (CS) techniques. An example of application of CS is shown in [18], in which different multiplexing alternatives to OFDM are analyzed and the CS technique is employed to overcome the random assignment of subcarriers. Indirect estimation methods, instead, require estimating the communication channel and, consequently, compensating for the contribution due to channel symbols (e.g., see [2, eq. (20)]), which are known at the RX side of any colocated radar. Indirect sensing methods can be divided in: 1) discrete Fourier transform (DFT)-based or correlation-based methods (i.e., methods based on the matched filter, MF, concept) [10], [12], [19]; 2) subspace methods [15], [16]; 3) maximum likelihood (ML)-based methods [14], [20], [21]. Correlation-based and DFT-based algorithms are conceptually simple and computationally efficient, but can generate poor radar images in the presence of closely spaced targets and/or strong clutter [22]. Moreover, they can be outperformed by subspace methods, like the wellknown multiple signal classification (MUSIC) algorithm and the estimation of signal parameters via rotational invariant technique (ESPRIT) at the price, however, of a significantly larger computational complexity [16]. An estimation accuracy comparable to that of subspace methods is provided by various ML-based algorithms, which also require a significant computational effort; relevant contributions to this field can be found in [14], [20] and [21].
Indirect sensing methods for JCAS MIMO systems employing OFDM are investigated in [4], [6], [17], [23], [24], [25], [26], [27], [28] and can be classified according to the same criteria as defined above for SISO systems. Correlation-based or DFT-based methods are developed in [4], [6], [17] and [28]. In particular, an algorithm combining matched filtering for range estimation with a onedimensional (1D)-MUSIC algorithm for both Doppler and azimuth (i.e., angle of arrival, AoA) estimation is proposed in [6]. In [4] a constant false alarm rate (CFAR) technique is employed to detect multiple targets, whereas the estimation of their range and Doppler is based on a two-dimensional (2D) fast Fourier transform (FFT); note that in this case the estimation of AoA is ignored and the availability of antenna arrays is beneficial for communication only. In [17] a 2D-FFT technique is exploited to estimate the range of multiple targets and their azimuth, whereas in [28] multiple FFTs are combined with a clutter removal technique to estimate their range, azimuth and Doppler.
Subspace methods are investigated in [23], [24] and [25]. More specifically, in [23] a subspace-based algorithm assisted by CFAR pre-processing is developed to estimate the range, velocity, and azimuth of multiple targets on the basis of a reduced number of samples and without resorting to high-resolution spectral estimation. In [24] the estimation accuracy of the 2D-MUSIC algorithm is assessed for a varying number of available data snapshots, whereas in [25] the use of an augmented beam-space approach is proposed to make the use of 2D-MUSIC and 2D-ESPRIT possible when hybrid digital arrays are used.
Methods based on an ML approach are proposed in [26] and [27]. In particular, the strategy devised in [26] is based on: 1) a preliminary channel estimation; 2) an ML-based technique for the estimation of the range and direction of arrival (DoA) of multiple targets. In [27], instead, a 1D technique, based on a systematic phase correction method and leveraging on the virtual array concept, is derived for the estimation of the azimuth of a single target. Finally, an hybrid approach, combining DFT-based and subspace methods, is illustrated in [29]. More specifically, channel estimation is accomplished, after DFT-based processing for range estimation, by means of an amplitude and phase estimation (APES)-based method; this method, that exploits the presence of inter-carrier interference (ICI) to produce a preliminary estimate of target Doppler, is followed by the 1D-MUSIC algorithm for angular estimation and by 1D-FFT processing to generate the final estimate of target Doppler.
It is important to point out that all the methods mentioned above for JCAS MIMO systems represent partial solutions to the problem of jointly estimating the range, Doppler, azimuth and elevation of multiple targets, i.e., briefly, to the four-dimensional (4D) imaging problem considered in this manuscript. On the one hand, the authors of [4] focus on range and Doppler estimation only; moreover, they develop a correlation-based method, which is employed after channel estimation and time-frequency synchronization in the context of long-range radar processing. On the other hand, the correlation-based or subspace methods proposed in [17], [24], [25] are able to compute accurate estimates of target azimuth only, whereas the other target parameters are ignored (see [24] and [25]) or estimated with limited accuracy, without any refinement process (see [17]). Similar considerations apply to the ML-based method illustrated in [27], since it can provide accurate estimates of the azimuth of multiple targets, but coarse estimates of their range only. These limitations originate from the fact that an acceptable complexity in subspace-based and ML-based methods can be achieved by neglecting (or accepting poor accuracy in) the estimation of a portion of target parameters. These considerations have motivated the work described in this manuscript, which aims at illustrating how a family of novel and computationally efficient sub-optimal methods able to generate accurate 4D radar images in an OFDM-based JCAS system equipped with TX and RX arrays can be developed; therefore, unlike previous research work, this paper aims at providing a full solution to the above mentioned 4D imaging problem. More specifically, the contribution provided by our manuscript is threefold and can be summarized as follows: 1) A general strategy, called Doppler-range-angle estimation with successive compensation (DRAEC), is proposed for the detection of multiple targets and the estimation of their parameters in a MIMO OFDM-based JCAS system. This strategy is based on the idea of: a) decoupling the problem of range-Doppler estimation from that of azimuth-elevation estimation; b) employing an algorithm for the estimation of the parameters of 2D complex tones to solve each of these two sub-problems.
2) An overview of the algorithms that can be employed for the detection and the estimation of 2D complex tones, and an analysis of their computational complexity are provided. On the one hand, five of these algorithms are known DFT-based estimation techniques, that have been proposed in the technical literature for related applications (namely, for harmonic retrieval [30], [31] or for radar sensing applications [32], [33]). On the other hand, the remaining algorithm, called extended Lee algorithm (ELA), is new, even if it can be considered as an extended version of the ML-based 1D algorithm developed in [27]. Note that all these algorithms have been originally proposed to solve estimation problems that are formally different from the ones considered in this manuscript. For this reason, in this manuscript, their use for target range-Doppler estimation and azimuth-elevation estimation according to the DRAEC strategy is illustrated in detail. Moreover, a unified notation is adopted in their description to ease their implementation and the reading of this manuscript.
3) Based on the above-mentioned estimation algorithms, seven different embodiments of the DRAEC strategy are proposed and compared in terms of both accuracy and computational complexity. Six of them are based on the known DFT-based estimation techniques mentioned above, whereas the remaining one relies on the CSFDEC algorithm. This allows us to assess how various state-of-the-art algorithms perform in 4D imaging and how large is the computational effort they require.
The remaining part of this manuscript is organized as follows. In Section II, the processing accomplished in a MIMO OFDM-based radar system is summarized and the received signal model adopted in our work is briefly derived. The DRAEC strategy is illustrated in Section III, whereas various estimators of 2D complex tones and their computational complexity are described in Sections IV and V, respectively. Different embodiments of the DRAEC strategy are proposed and compared, in terms of accuracy and complexity, in Section VI. Finally, some conclusions are offered in Section VII.
Notation: Throughout this paper, the following notation is adopted: 1) (·) T denotes matrix transposition; 2) (·) * and (·) H denote complex conjugate and complex conjugate transpose (Hermitian operator), respectively; 3) the symbols ⊗, , * and × represent the Kronecker, Hadamard, Khatri-Rao and Cartesian product operators, respectively; 4) {x} and {x} indicate the real part and the imaginary part, respectively, of the complex variable x; 5) diag(v) represents a square diagonal matrix having the elements of the vector v along its main diagonal; 6) 0 M,N denotes the M × N null matrix; 7) I M denotes the order M identity matrix; 8) Ø represents the empty set.

II. SYSTEM AND SIGNAL MODELS
This section focuses on the architecture of the MIMO OFDM-based JCAS system considered in our manuscript and on the processing accomplished at its receive side. Our main objectives are deriving the received signal model in the presence of multiple targets and illustrating some essential assumptions on which it relies. The architecture of the considered JCAS system is illustrated in Fig. 1 and has the following essential features: 1) Its transmitter is colocated with the receiver; consequently, the receiver has full knowledge of the structure and content of the transmitted signal and of its carrier frequency, and exploits these information for sensing purposes only.
2) It is equipped with a transmit (TX) horizontal uniform linear array (HULA) and a receive (RX) vertical uniform linear array (VULA), consisting of N T and N R elements, respectively. All the antennas are placed on the same planar shield, so that a 2D reference system lying on the plane of the physical antenna array can be defined, as illustrated in Fig. 2 (where λ denotes the wavelength of the transmitted signal).
3) The data frames it transmits are made of M consecutive OFDM symbols, each consisting of N subcarriers. Such symbols can convey both pilot tones (for channel estimation and synchronization) and information data to be sent to a single or multiple receivers at different locations.
In our work, each couple of physical TX and RX antennas is replaced by a single virtual antenna (VA); the abscissa x (p,q) v and the ordinate y (p,q) v of the VA element associated with the pth TX antenna and the qth RX antenna (briefly, the (p, q) VA) are evaluated as (e.g., see [34, eqs. (1)-(2)]) and respectively, with p = 0, 1, . . . , N T − 1 and q = 0, 1, . . . , r ) denote the coordinates of the pth TX and qth RX antenna, respectively. It is easy to show that the set of N VA = N T N R VAs associated with the physical arrays shown in Fig. 2 forms a virtual uniform rectangular array (URA). Moreover, based on (1) and (2), the abscissa and ordinate of the VA (p, q) are and r ) represent the coordinates of the leftmost TX (lowermost RX) antenna and d t (d r ) denotes the distance between adjacent antennas of the TX (RX) array (see Fig. 2, where it is assumed that d t = d r = λ/2).  In the following derivations, we concentrate on the transmission of a single frame. The complex envelope of the radio frequency (RF) signal conveying the mth OFDM symbol radiated by the pth TX antenna (with p = 0, 1, . . . , N T − 1) can be expressed as (e.g., see [14, eq. (3)]) up to a transmit delay; here, q(t) is a windowing function, s m,n is the mth channel symbol conveyed by the nth subcarrier and transmitted by the pth TX antenna, f = 1/T is the subcarrier spacing, T is the OFDM symbol interval, T s T + T G is the overall duration of the OFDM symbol and T G is the cyclic prefix interval. Following [14], a rectangular windowing function is adopted, so that q(t) = 1 for t ∈ [−T G , T] and q(t) = 0 elsewhere. Given the complex envelope (5), the RF waveform radiated by the pth TX antenna in the considered frame can be expressed as where f c = c/λ denotes the frequency of the local oscillator employed in the up-conversion at the TX side and c the speed of light. Let assume now that x (p) (6) is reflected by K distinct point targets, and that the kth target (with k = 0, 1, · · · , K−1) is located at the (initial) distance R k from the transmitter, moves with the radial velocity 1 v k with respect to it and is characterized by the azimuth (elevation) angle θ k (φ k ). It is not difficult to show that the complex envelope of the RF signal r (q) Fig. 1) captured by the qth RX antenna is 2 (e.g., see [14, eqs. (5) and (6)]) 1. This velocity is positive (negative) if the target approaches (moves away from) the JCAS system.
2. Note that the overall delay that characterizes the echo originating from the kth target depends on all the parameters of the target itself (namely, its range, its velocity and its angular parameters) and changes over time.
where (8) and α (p,q) k are the overall propagation delay and the attenuation, respectively, associated with the kth point target (and observed on the VA (p, q)), is the Doppler shift due to the motion of the kth target and w (q) (t) is the complex additive white Gaussian noise (AWGN) process affecting r (q) (t).
The signal r (q) (t) in (7) undergoes analog-to-digital conversion followed by DFT processing. A simple mathematical model that describes the sequence generated by sampling r (q) (t) in the mth OFDM symbol interval can be derived as follows. Substituting the right-hand side (RHS) of (5) in that of (7), extracting the portion associated with the mth OFDM symbol from the resulting expression and substituting t with t = t − mT s yields where and Note that: 1) the term D(f , t ), in (10), is responsible for the so called range migration effect (see [35]) due to the kth target Doppler frequency, whereas ν m (f D k ) is proportional to both f D k and the OFDM symbol index m; 2) the phase of A(R k ) in (10) depends on the range R k only, whereas that of γ n (R k ) in (14) is proportional to both R k and the subcarrier index n (see (13) and (16), respectively); 3) The terms in (15) and in (17) both depend on the virtual array configuration as well as the target DoA; 4) the term ξ n (f D k , t ) in (14) produces a time-dependent phase rotation influenced by both the Doppler frequency f D k and the subcarrier index n (see (18)); 5) the term ζ m,n (f D k ) in (14) introduces a phase rotation depending on both the OFDM symbol index m and the subcarrier index n (see (19)), and accounts for the so-called intersubcarrier Doppler effect (e.g., see [14, Sec. II, p. 3]). It is not difficult to show that sampling r (q) m (t ) (10) at the instant t l lT/N (i.e., sampling r (q) (t) at the instant t l t l + mT s ), with l = 0, 1, . . . , N − 1, yields 3 here, the term D l (f ) exp(j2π flT/N) accounts for the ICI effect due to the range migration, is the Gaussian noise affecting r  If the target Dopplers {f D k } are sufficiently small and, more precisely, |f D k /f c | 1/(MN) for any k, the factors ξ n,l (f D k ) and ζ m,n (f D k ) appearing in the RHS of (21) can be neglected; this leads to the simplified signal model The N signal samples acquired in the mth OFDM symbol interval through the qth RX antenna undergo serial-toparallel (S/P) conversion, as shown in Fig. 1; this produces the N-dimensional vector r where W (q) n is the AWGN sample affecting the nth subcarrier. The received signal model expressed by (24) is general, but quite complicated. In our work, a simplified version of it can be employed since: 1) the approximation is made for any n, since f c is assumed to be much greater than f (see (17)); 2) it is assumed that, in the transmission of each OFDM symbol, disjoint subsets of the N available subcarriers are assigned to distinct TX antennas.
The last assumption means that the signal radiated by each TX antenna can be represented as the superposition of a set of multiple subcarriers exclusively assigned to that antenna, i.e., briefly, it consists of private subcarriers only. 4 Therefore, if S (p) m denotes the set collecting the indices of the subcarriers radiated by the pth TX antenna in the mth symbol interval, we have that for any p 1 = p 2 , with p 1 and p 2 ∈ {0, 1, . . . , N T − 1}, and for any m, where J {0, 1, . . . , N − 1}. It is important to point out that, although our assumption about the use of subcarrier frequencies is strong and entails a reduction by a factor N T in the transmission rate (with respect to the case in which all the subcarriers are employed by each TX antenna), it paves the way for the development of target detection and estimation algorithms requiring a limited computational effort. In fact, under the last assumption, the sum over p appearing in the RHS of (24) involves a single term different from zero; consequently, (24) becomes which represents a sample of the 2D signal employed for channel estimation. Here, p a denotes the index of the TX antenna to which the nth subcarrier has been assigned in the mth OFDM symbol interval (the dependence of p a on 4. The use of private subcarriers in JCAS systems has been first proposed in [19] to improve the accuracy of target range and angle estimation through a CS technique. m and n is not explicitly shown to ease notation). Based on (3), (4), (11)-(13), (15) and (16), it is easy to show that (30) and where with z = p a , q, m or n, X = H k , V k , D k or ρ k , and denote the normalized horizontal frequency and the normalized vertical frequency, respectively, associated with the kth target, is the normalized Doppler frequency, is a normalized frequency accounting for both the Doppler of the kth target and its range through the normalized target delay 5 and Note that: 1) F ρ k and F D k satisfy the inequalities F ρ,min ≤ F ρ k ≤ F ρ,max, and F D,min ≤ F D k ≤ F D,max , respectively, with F ρ,min = 0, F ρ,max = 1, F D,min = −1/2 and F D,max = 1/2 for any k; 2) F V k and F H k satisfy the inequalities −d r /λ ≤ F V k ≤ d r /λ and −d r cos(φ k )/λ ≤ F H k ≤ d r cos(φ k )/λ, respectively, for any k; 3) the range of F H k is also limited by the elevation φ k of the kth target for any k (see (33)); 4) the ranges of F V k and F H k are maximized for d r = λ/2 and d t = λ/2, respectively.
Based on (29)-(31), eq. (28) can be rewritten as 6 5. Note that, F r k is always positive, whereas F D k is positive (negative) if the kth target is approaching (moving away from) the considered JCAS system. 6. Note that channel symbols {s denotes an estimate of the channel frequency response H (p a ,q) m,n characterizing the nth subcarrier frequency in the mth OFDM symbol interval for the VA (p a , q), A k α k exp(jω k ) and is the noise sample affectingĤ m,n , in (40), represents the noisy measurement produced by the symbol compensation block for the VA (p a , q), the nth subcarrier and the mth OFDM symbol interval (see Fig. 1, where this block is followed by a buffer storing the measurements acquired over each frame). Moreover, the measurements acquired over the OFDM frame and the whole virtual array are collected in the set 8 m,n ] represents the M × N matrix formed by all the measurements acquired through the VA (p, q) over the OFDM frame. Note that this set consists of N VA matrices, one for each of the N VA VAs, and that it represents the output of the last block appearing in Fig. 1.
The measurement model (40) deserves the following comments: 1) The complex gain A k appearing in its RHS accounts for the phase rotation due to the path delay, the path loss and the gain (attenuation) introduced by the kth target.
2) Since, in all our computer simulations, a N s ary phase shift keying (PSK) constellation is adopted for the channel symbols {s n } (see (28)), i.e., it can be modelled as AWGN (the variance of each sample is denoted σ 2 W ).

3) The noisy samples {Ĥ
m,n } of the 4D channel response acquired over a single OFDM frame can be modelled as the superposition of an AWGN process with K 4D complex exponentials, whose amplitude, phase and frequencies provide information about the range, the Doppler and the angular coordinates of the detectable targets. For this reason, target detection and estimation is tantamount to identifying the K complex exponentials that form the useful component 7. In the following, the dependence of α (p a ,q) k on (p a , q) is neglected; therefore, α k = α (p a ,q) k is assumed. 8. In the following, the antenna index p a is replaced by p (with p = 0, 1, . . . , N T − 1) for simplicity. of the sequence {Ĥ (p,q) m,n } and to estimating their parameters, respectively.
Finally, it is important to make some considerations about the use of private subcarriers and the criteria that can be adopted in the selection of their subsets {S (p) m }. From the rules illustrated above about the use of private subcarriers it can be easily inferred that the maximum data rate achievable through the proposed transmission scheme is identical to that of a JCAS system equipped with a single TX antenna, i.e., as already mentioned above, it is N T times lower than that provided by a MIMO system with shared subcarriers (e.g., see [19]). In addition, the TX array is not exploited for beamsteering, as suggested, for instance, in [29], where a single-stream beamforming model is assumed. Despite this, the considered JCAS system benefits from the availability of a TX array, since this results in a larger virtual array (i.e., in an increase of the overall number of VAs, N VA ) and, consequently, in a better angular resolution [27]. As far as the selection of the subsets {S (p) m } is concerned, in our computer simulations, a pseudo-random mechanism has been adopted in assigning the N available subcarriers to the N T TX antennas in the transmission of the mth OFDM symbol of a given frame; moreover, the same pseudo-random pattern has been employed for all the transmitted frames. The choice of this strategy is motivated by the fact that randomly changing the subset of subcarriers from symbol to symbol allows the considered radar system to benefit from transmit diversity.

III. DESCRIPTION OF THE PROPOSED APPROACH TO THE ESTIMATION OF MULTIPLE TARGETS
In this section, the problem of developing reduced complexity methods for the detection of multiple targets and for the estimation of their parameters in the MIMO OFDM-based JCAS system described in the previous section is tackled. We first describe a general strategy to devise novel solutions to this problem. Then, we provide some indications about the processing to be accomplished by each of the two main parts it consists of.

A. DESCRIPTION OF THE PROPOSED STRATEGY
Achieving joint ML estimation of an unknown number of targets, given the set of measurements S H , in (42), is an overly complicated problem, since it involves a large number of parameters to be estimated (more precisely, five parameters per target plus the overall number of targets), even for small values of K. This motivates our interest in the development of sub-optimal methods based on the idea of turning a multidimensional estimation problem into a set of interconnected lower dimensional sub-problems. In the remaining part of this subsection, we illustrate a general strategy, called DRAEC, for the derivation of a new class of such methods. According to this strategy, range and Doppler estimation is decoupled from angular estimation. This explains why its structure, described by the block diagram shown in Fig. 3, contains two core blocks, called range-Doppler estimator (RDE) and angular estimator (AE): in fact, the former block accomplishes target detection and jointly estimates target range and Doppler, whereas the second one identifies multiple targets characterized by similar Doppler and ranges, and estimates their azimuth and elevation. Note also that the proposed structure includes a fusion block, whose task is merging the information provided by the first two blocks in order to generate a 4D radar image in the form of a point cloud. The processing accomplished by the core blocks can be summarized as follows. Based on the available measurements (i.e., on the set S H , in (42)), the RDE generates the so-called target range-Doppler profile (TRDP), namely a collection of: 1) range-Doppler couples at which relevant echoes are detected; 2) an estimate of the complex amplitude associated with each of these couples (the absolute value of such amplitudes allows us to rank the couples on the basis of their perceptual importance). More precisely, the TRDP is represented by the set whereF D k ,F r k andÂ k denote the estimates of the normalized Doppler frequency, the normalized delay and the complex amplitude, respectively, associated with the kth range-Doppler bin 9 in which (at least) one target has been detected, 10 and K (RDE) is the overall number of relevant range-Doppler bins identified by the RDE. The set S (RDE) , in (43), is passed to the AE, which processes it jointly with the set S H (see (42)) in order to: 1) identify all the targets associated with each of the range-Doppler-complex amplitude triplets forming the TRDP; 2) generate the so-called angular profile (AP), that collects the estimates of the angular parameters and the complex amplitude of all the targets detected within each range-Doppler bin. In practice, the AP information associated with the kth element (i.e., triplet) of S (RDE) is represented by the set denote the estimates of the normalized horizontal frequency, normalized vertical frequency and complex amplitude of the lth target detected in the kth range-Doppler bin, respectively, whereas K (AE) k represents the overall number of targets detected in that bin. The normalized frequenciesF H k [l] and F V k [l] are jointly processed by the fusion block in order to generate, on the basis of (33) and (34), the estimateŝ θ k [l] andφ k [l] of the azimuth θ k [l] and the elevation φ k [l], respectively, characterizing the lth target identified in the kth range-Doppler bin. Moreover, the fusion block processes the 9. As shown in the next section, the 2D FFT processing executed by the RDE leads to discretizing the range-Doppler domain and, in particular, to partitioning it into multiple range-Doppler bins.
10. Note that the RDE is unable to separate multiple targets whose range and Doppler fall in the same bin. estimatesF D k andF r k to compute the estimatesv k [l] and R k [l] of the velocity v k [l] and range R k [l], respectively, on the basis of (9), (35) and (37). Then, the above mentioned estimates are collected in the set where represents the final output of the DRAEC strategy if the RDE is not exploited again. Alternatively, it can be passed to the RDE with the aim of re-estimating the range and Doppler parameters of each of the K (AE) k detected targets for any k; this step is expected to generate an updated version of the set S (RDE) (see (43)); fusing this set with the overall AP produces an updated version of the set S (DRAEC) , that collects finer estimates of the parameters of all the detected targets.

B. PROCESSING TASKS ACCOMPLISHED BY THE CONSTITUENT BLOCKS
In this subsection, the essential processing tasks accomplished by the RDE and AE blocks appearing in Fig. 3 are sketched; in our description it is assumed, without any loss of generality, that the considered MIMO OFDM-based JCAS system is equipped with the URA illustrated in Fig. 2.
The RDE extracts from the set S H (see (42)), collecting N VA matrices, a single matrix, denotedĤ (p R ,q R ) and referring to a specific VA (called reference VA and associated with the choice (p, q) = (p R , q R ); see Fig. 2). Then, it processeŝ H (p R ,q R ) to generate the TRDP, i.e., the set S (RDE) , in (43). It is important to point out that: 1) Based onĤ (p R ,q R ) , the TRDP can be generated by estimating the complex amplitudes and the frequencies of the complex exponentials that form the useful component of the 2D sequence {Ĥ 2) The parameter K (RDE) appearing in (43) represents the overall number of range-Doppler bins in which at least a single target is detected. In fact, at this stage, the RDE detects multiple targets, characterized by similar ranges and Dopplers (and, in particular, such that their parameters fall inside in the same range-Doppler bin) as a single target. This explains why, in general, the kth range-Doppler bin selected by the RDE (with k = 0, 1, . . . , K (RDE) − 1) may contain multiple (say, K (AE) k ) targets, which are overlapped in the range-Doppler domain, but have distinct angular coordinates.
3) The absolute value of the complex amplitudeÂ k appearing in (43) represents the perceptual importance of the kth detected range-Doppler bin (with k = 0, 1, . . . , K (RDE) ); in the following, we assume that the elements of S (RDE) are ordered according to a decreasing perceptual importance, so The set S H and the TRDP (see (42) and (43), respectively) feed the AE, that sequentially accomplishes the three steps listed below.
1) It checks if the kth range-Doppler bin satisfies the inequality with k = 0, 1, . . . , K (RDE) ; here, T (RDE) is a proper threshold. Any bin not meeting this condition is discarded. This leads to the reduced TRDP 2) It merges the information provided by the N VA matrices of the set S H (see (42)) inK (RDE) N T × N R matrices, one for each ofK (RDE) range-Doppler bins selected in the previous step. The kth matrix (with k = 0, 1, . . . , ]; the element appearing on its pth row and qth column is evaluated as with p = 0, 1, . . . , N T − 1 and with q = 0, 1, . . . , N R − 1; here,F ρ k represents the estimate of the normalized frequency in (36).
3) It processes the matrixH k to identify all the targets contained in the kth range-Doppler bin in order to generate the set S The AE output, i.e., the overall AP, results from merging all the information contained in theK (RDE) sets {S (AE) k } and collects the estimates of all the angular parameters referring to distinct targets. This concludes the AE processing.
It is worth noting that: ) appearing in the RHS of (49) aims at compensating for the factor a m (F D k ) (a n (−F ρ k )) which is visible in the same side of (40); in other words, it is expected to cancel the dependence ofH (p,q) k on the Doppler and range of the kth target (leaving, however, the dependence on its normalized horizontal and vertical frequencies).
2) Similarly as step 1) of the RDE, step 3) of the AE processing requires estimating the complex amplitudes and the frequencies of the overlapped complex exponentials that form the useful component of the 2D sequence {H If the AE output is passed to the RDE, the last block sequentially accomplishes the three steps 11 illustrated below for each of the consideredK (RDE) range-Doppler bins; in our description of such steps we refer to the kth range-Doppler bin (with k = 0, 1, . . . ,K (RDE) − 1).
1) The RDE generates the subset where T (AE) is a proper threshold.
2) It merges the information provided by the N VA antennas at the N subcarrier frequencies and over the M OFDM symbol intervals (i.e., over the whole OFDM frame) to generate a set ofK (AE) k M × N matrices, each referring to a single element of the setS (AE) k . More specifically, the matrix associated with the lth element ofS [Ȟ m,n [k, l]]; moreover, the element appearing on its mth row and nth column is evaluated aš ) of the normalized Doppler frequency, normalized delay and complex amplitude, respectively, of the above mentioned target. These information are collected in the set 11. Actually, the first step is ignored if all the targets detected in each range-Doppler bin are taken into consideration.
Note that the term a p (F H k [l]) (a q (F V k [l])) appearing in the RHS of (53) aims at compensating for the factor a p (−F H k ) (a q (−F V k )) that appears in the same side of (40); in other words, it is expected to cancel the dependence ofĤ m,n [k, l] on the normalized horizontal and vertical frequencies (i.e., on the angular parameters) of lth target detected in kth range-Doppler bin.
The strategy described above is called Doppler-rangeangle estimation with successive compensation, or DRAEC, and is summarized in Algorithm 1. Its final output is represented by the set that collects the estimates referring to represents the contribution due to all the targets detected in the kth range-Doppler bin. Note that, the setS (DRAEC) , in (55), results from merging all the information provided the RDE and the AE, namely the setsS (RDE) and {S (AE) k }. In other words, the final output is obtained by 1) converting each normalized frequency to the corresponding spatial parameter; 2) associating each target with its set of spatial parameters. These operations are carried out by the fusion block appearing in Fig. 3.
Finally, it is important to stress that the most important task accomplished by both the RDE and the AE is represented by the estimation of the parameters of 2D complex oscillations on the basis of a set of noisy measurements. Different algorithms can be employed for this task; the selection of a specific algorithm leads to a different instance of the DRAEC strategy, as illustrated in the following sections.

IV. DESCRIPTION OF VARIOUS ALGORITHMS FOR THE DETECTION AND THE ESTIMATION OF TWO-DIMENSIONAL COMPLEX OSCILLATIONS
In this section, we concentrate on the problem of detecting multiple overlapped 2D complex exponentials and estimating their parameters in the presence of AWGN. Various estimators, representing distinct solutions to this problem, are illustrated. In all cases, essential mathematical details are provided and a unified mathematical notation is employed. More specifically, it is assumed that: represents the input sequence 12 of any 2D estimator employed in the RDE or in the AE; 2) the element (m, n) of this sequence can be expressed aŝ where have the same meaning as the corresponding parameters appearing in the RHS of (40) and W (X) m,n is the noise sample affectingĤ (X) m,n (an AWGN model is adopted for the sequence {W (X) m,n }). All the algorithms described in the following subsections make use of 2D periodograms. More specifically, target detection and range & Doppler estimation in the RDE require the computation of the M 0 × N 0 matrix (59), is given by where and L (RDE) 1 and L (RDE) 2 represent the oversampling factors adopted in RDE processing. Note that Y (RDE) [l, p] is associated with the normalized Doppler frequency and the normalized frequency (accounting for both range and Doppler; see (36)) Similarly, angular estimation in the AE requires the computation of theM 0 ×N 0 matrix and L In the remaining part of this section we take into consideration six different estimation algorithms and provide a brief mathematical description of each of them. The first five algorithms are FFT-based methods and, more precisely, are the 2D periodogram method [2], the CSFDEC algorithm derived in [32], the estimation algorithm proposed by Popović-Bugarin and Djukanović in [30], a modified version of the estimation algorithm devised by Fan et al. in [33] for solving the problem of channel estimation in a hybrid millimeter-wave massive MIMO system and the q-shift estimator (dubbed QSE) developed in [31]. The sixth (and last) algorithm can be considered as an extension to 2D frequency estimation of the algorithm developed by Lee and Chun in [27]; therefore, it is dubbed extended Lee algorithm, or ELA. In the description of each estimation algorithm, we first describe its formulation for the RDE; then, we illustrate the changes required to make its use possible in the AE.

A. TWO-DIMENSIONAL PERIODOGRAM METHOD
This method is based on the idea that the frequencies of the 2D complex exponentials forming the useful component of the sequence {Ĥ (X) m,n } are associated with the peaks of a 2D periodogram. For this reason, if the kth target is considered (with k = 0, 1, . . . ,K−1, whereK denotes an estimate of the overall number of targets), the estimatesF D k andF ρ k of its normalized frequencies F D k (35) and F ρ k (36), respectively, are evaluated asF (64) and (65), respectively), where (l ) is the value of the couple (l,p) associated with the kth of thê K most relevant local maxima (peaks) of the 2D sequence for any positive integer U. This algorithm can be also used in the AE block to evaluate the estimatesF (72) and (73), respectively) of the normalized frequencies F H k (33) and F V k (34); the couple of (l The estimation accuracy of this method can be improved by: 1) Extracting an I l × I p sub-matrix (where I l and I p denote the interpolation orders adopted in the Doppler and range domains, respectively), whose central element is k ] (with k = 0, 1, . . . ,K − 1 and X = RDE or AE) from Y (X) in order to generate a more detailed representation of the analyzed spectrum through the interpolation 13 of the sub-matrix elements.
2) Identifying the peak of the interpolated spectrum over the considered 2D domain.

B. COMPLEX SINGLE FREQUENCY DELAY ESTIMATION AND CANCELLATION ALGORITHM
The second FTT-based method combines a single 2D tone estimator, named complex single frequency delay estimator (CSFDE), with a serial cancellation procedure. This algorithm can be exploited by the RDE block as it is, whereas some modifications are required if it is employed in the AE block. For this reason, we first focus on its use in the RDE and, then, we illustrate the changes to be made for its use in the AE. 13. In our computer simulations, the 'spline' interpolation (interp2 function) of MATLAB R2022b has been employed.
The derivation of the CSFDE is based on the assumptions that: 1) K = 1 in the model (58) forĤ (RDE) m,n ; 2) the unknown frequencies F 1,0 = F D and F 2,0 = −F ρ appearing in that model can be represented as F D = F D,c + δ DFD and F ρ = F ρ,c + δ ρFρ , respectively, where F D,c (F ρ,c ) is a coarse estimate of F D (F ρ ), δ D (δ ρ ) is the associated residual, andF D (F ρ ) is expressed by (66) (67) and represents the normalized fundamental frequency characterizing an order M 0 FFT (N 0 IFFT).
The CSFDE first exploits the 2D periodogram method to detect a single tone and to compute a coarse estimate of its complex amplitude A and frequencies; then, it makes use of an iterative procedure for estimating the residuals and refining the complex amplitude. More specifically, in the initialization of the CSFDE, the following quantities are computed: 1) The set of is the order (M 0 , N 0 ) DSFT (see (59)) of the matrixĤ (k 1 ,k 2 ) ZP that results from zero padding 14 with m = 0, 1, . . . , M − 1 and n = 0, 1, . . . , N − 1.
2) The coarse estimatesF 3) The initial estimatê of the complex amplitude A; here, 15 14. In practice, the structure of the matrixĤ and under the assumption that with X = and X = , respectively. 6) The initial fine estimateŝ of F D and F ρ , respectively. This concludes the initialization, which is followed by an iterative procedure whose index i is set to 1. The ith iteration of this procedure is fed by the estimatesF are computed. It is worth mentioning that an alternative to (81) for the evaluation ofȲ k 1 ,k 2 (F D ,F ρ ) in (82)-(85) is represented by the use of a 2D interpolation method applied to a I l ×I p submatrix 16 ofȲ k 1 ,k 2 (77); the need of interpolation originates 16. In the following, I M and I N denote the interpolation orders adopted for the first and second dimension, respectively, of that matrix. 2) Estimation of the complex amplitude -The new estimatê A (i) ofÂ is evaluated by means of (80); in doing so, the couple (F (i) D ,F (i) ρ ) is used in place of (F (0) D,c ,F (0) ρ,c ). After that the last step has been carried out, the index i is incremented by one and a new iteration is started. At the end of the last (i.e., of the N it th) iteration, the fine estimateŝ respectively, become available and the algorithm stops.
The CSFDE algorithm represents the core of the CSFDEC algorithm, which is used to recursively estimate the multiple tones forming the useful component of the complex sequence {Ĥ (RDE) m,n }, whose (m, n)th element is expressed by (58) with K ≥ 1 and, in general, unknown. The CSFDEC algorithm is initialized by: 1) Running the CSFDE algorithm to compute the initial of the parameters F D 0 , F ρ 0 and A 0 , respectively, that characterize the first target.
2) Setting the recursion index i to 1 andȲ Then, a recursive procedure is started. The ith recursion of this procedure is fed by the vectorsF  2) Multiple (say,N it ) iterations are executed to refine the estimate of the parameters of the new tone detected 17. If K is known, the computation of the energy ε 0,0 [i] and its comparison with a threshold are not required; in fact, in this case, the CSFDEC algorithm stops at the end of its Kth recursion.
in the previous step. The processing accomplished in this step follows closely that described in the refinement part (i.e., in the second step) of the CSFDE. For this reason, in each iteration, a new estimate of the complex amplitude and of the two residuals of the ith tone are computed.
3) Each of the i detected tones is re-estimated after cancelling the leakage due to all the other (i − 1) tones. This allows to progressively refine the amplitude, normalized Doppler frequency and normalized delay of each tone, thus generating the final estimates. Note that, in principle, this re-estimation procedure can be repeated multiple (say, N REF ) times.
As already mentioned above, the CSFDEC algorithm can also be employed in the AE. However, this requires: 1) UsingĤ 2) Evaluating the initial estimate of the complex amplitude of the strongest 2D tone as (see (80)) (72) and (73), respectively) are the coarse estimates of the normalized vertical and horizontal frequencies, respectively.
3) Replacing (84) and (85) and respectively. The CSFDEC algorithm is summarized in Algorithm 2 for the case in which K is unknown. If K is known, in step d) of the frequency estimation procedure, the evaluation of the energy ε 0,0 [i] and its comparison with a threshold are replaced by a comparison of the iteration index i with K; if i = K, step e) is executed, otherwise i is increased by one and the algorithm proceeds with step b).

C. POPOVIĆ ALGORITHM
The third FFT-based method (namely, Alg-P) computes the frequency estimates through a serial refinement and cancellation procedure based on the computation of a set of shifted DFT coefficients and their subsequent parabolic fitting. If employed in the RDE, it is fed by the complex sequence {Ĥ  (0, 3), (3, 0) and (3, 3). Then, set the recursion index i to one and the re-estimation index r to one. 2 Frequency estimation: b-Run the CSFDE algorithm to compute the estimateŝ F (i) and n = 0, 1, . . . , N − 1. Then, it sequentially executes the three steps described below for the kth target (with k = 0, 1, . . . ,K − 1, whereK is an estimate of K).
1) 2D periodogram maximization -In this step, the coarse estimates of the normalized Doppler frequency F D k and the normalized range frequency F ρ k are evaluated asF k ], respectively (see (64) and (65), respectively); here, m,n is illustrated below).
2) Frequency refinement and amplitude estimation -In this step, the fine estimatesF D k andF ρ k of F D k and F ρ k , respectively, are computed according to the formula (see [30, with X = D or ρ; here, NUM = ψ 2 X,3 P X,1 − P X,2 + ψ 2 X,2 P X,3 − P X,1 + ψ 2 X,1 P X,2 − P X,3 , (96)  Table 2, eq. (12)-(13)]). If the energy ε[k + 1] Y (k+1) 2 is smaller than T (P) , where T (P) is a proper threshold, the algorithm stops and the estimateK = k + 1 of K is generated; otherwise, k is increased by one and the three steps described above are accomplished again.
The use of Alg-P in the AE requires the following modifications: 1) the spectral coefficient Y (k) [l,p] appearing in the RHS of (94) is still expressed by (69), where, however,Ĥ (AE) m,n is replaced by 2) formula (95) (with X = V or H) is employed to compute the estimatesF H k andF V k of F H k and F V k , but the quantities and are used in place of P ρ,i (103) and P D,i (102), respectively. A schematic description of Alg-P is provided in [30, Sec. 3, Tables 1 and 2].

D. MODIFIED FAN ALGORITHM
The fourth FFT-based method (namely, the MFA) results from: 1) adapting the estimation algorithm devised in [33] to the signal model expressed by (58); 2) including zeropadding in the initialization of the 2D periodogram method (see Section IV-A). The processing accomplished by the proposed algorithm evolves through the following two consecutive steps; note that, in this case, an estimate, denoted K, of the overall number of targets is required. 1) 2D periodogram maximization -In this step the 2D periodogram is maximized (see our description of the 2D periodogram method in Section IV-A) in order to evaluate the coarse estimatesF D k ,c andF ρ k ,c of the normalized Doppler frequency F D k , in (35), and the normalized range frequency F ρ k , in (36), respectively, with k = 0, 1, . . . ,K − 1.
2) Frequency refinement -First, the fine estimates of normalized Doppler frequency F D k and the normalized range frequency F ρ k are evaluated aŝ with k = 0, 1, . . . ,K − 1,˜ D k (˜ ρ k ) represents the trial variable for the residual D k ( ρ k ), is the cost function selected for the considered estimation problem (see [33,Sec. III,eq. (31)]), and 1) The 2D periodogram method processes the spectral matrix Y (AE) , in (68), (in place of Y (RDE) , in (59)) to produce the coarse estimatesF H k ,c andF V k ,c of the normalized frequencies F H k (33) and F V k (34), respectively (with k = 0, 1, . . . ,K − 1).
2) The fine estimates of F H k and F V k are evaluated aŝ The estimates of the residualsˆ H k andˆ V k represent the solution of an optimization problem formally identical to (110), where, however, J k (˜ D k ,˜ ρ k ) is replaced by where of the complex amplitude characterizing the kth target is evaluated (with k = 0, 1, . . . ,K − 1).
The MFA is summarized in Algorithm 3.

E. Q-SHIFT ESTIMATOR
The fifth FFT-based algorithm (namely, the QSE) has been proposed by [31] to estimate the frequency of a 2D complex tone in the presence of AWGN. Similarly to the MFA, this estimator makes use of the 2D periodogram method for coarse frequency estimation and requires prior knowledge of the overall number of targets; however, it exploits a different method for frequency refinement. In fact, the last task is accomplished by a serial procedure that requires the evaluation of the DFT coefficients located at the relevant frequency bins shifted by a quantity q ∈ [−0.5, 0.5]. In practice, the final estimates of the normalized frequencies F D k (35) and F ρ k (36) are evaluated aŝ  ) is provided by the 2D periodogram method for the RDE (see Section IV-A). The estimation of the residuals (δ D k ,δ ρ k ), instead, is accomplished by an iterative procedure; this is initialized by setting the initial estimates of the residuals (namely,δ ρ k ) to zero and the iteration index i to 1. In the ith VOLUME 4, 2023 iteration (with i = 1, 2, . . . , N it , where N it is the overall number of iterations), the new estimate of the residualδ (i) with X = D or ρ, and k = 0, 1, . . . ,K − 1; here,K denotes our estimate of K, and are the DFT coefficients evaluated with small shifts (quantified by the real parameters q D and q ρ ) with respect to the periodogram peak associated with the couple and are the normalized frequencies associated with the shifts q D and q ρ , respectively, z is an integer belonging to the set {0, ±1}, and and are correction factors. The final estimates of the residuals (δ D k andδ ρ k ) are evaluated asδ D k =δ The use of the QSE in the AE requires the following modifications: 1) The final estimates of the normalized frequencies F H k and 2) The computation of the residualsδ H k andδ V k is still based on (120) (with X = H or V), but (122) and (121) are replaced by and respectively; here, are the normalized frequencies associated with the shifts q H and q V , respectively, and z is an integer belonging to the set {0, ±1}. Moreover, the correction factors c X (q X ) (with X = H or V) are both computed according to (126). A schematic description of the QSE is provided in [31, Sec. IV, Algorithm 2].

F. EXTENDED LEE ALGORITHM
The last algorithm (namely, the ELA) has been originally proposed in [27, Sec. III] to perform azimuth estimation in a MIMO radar equipped with a uniform linear array (ULA) and is based on an ML approach. However, in our work, the following modifications have been made: 1) The original algorithm, being developed to estimate the frequencies of 1D tones, has been adapted to the 2D signal model expressed by (58), so making its use possible in both the RDE and the AE.
2) An iterative procedure for frequency refinement has been added. In each iteration of this procedure, the grid adopted in the search for the frequency estimate of a given target is adjusted to improve the achieved accuracy.
3) The 2D periodogram method has been employed for the initialization of the ELA; in the original algorithm, instead, the frequency estimates are initialized to zero for all the detected targets.
The ELA is fed by the 2D sequence {Ĥ (RDE) m,n }; in its initialization, this sequence is processed by the 2D periodogram method to evaluate the initial estimatesF ρ k andÂ (0) k of F D k , F ρ k and A k , respectively, with k = 0, 1, . . . ,K − 1 (beingK a preliminary estimate of the overall number of targets K), and the iteration index i is set to one. Then, the refinement procedure is started. In its ith iteration (with i = 1, 2, . . . , N it , where N it is the overall number of iterations), the new estimatesF ρ k of F D k and F ρ k , respectively, are evaluated as with k = 0, 1, . . . ,K − 1; here, is the ML cost function evaluated for the trial couple This grid has the following relevant properties: 1) its center depends on bothF ; 2) its step sizes get smaller as i increases. More precisely, its node (z D , z ρ ) (with z D = 0, 1, . . . , N F D − 1 and z ρ = 0, 1, . . . , if F X,min ≤F if F X,min + δ X ≤F (i−1) if F X,max − δ X <F , so that, when i = 1, two adjacent bins of the spectrum considered in coarse frequency estimation are covered. 18. The dependence of I on the target index k is not shown in the following three equations to ease notation.

Algorithm 4: Extended Lee Algorithm (ELA)
Input: The matrixĤ (RDE) (Ĥ (AE) ) for the RDE (AE), the overall number of iterations in frequency refinement (N it ) and an estimate of the overall number of targets (K). 1 Initialization: Evaluate the estimates ( k ) by findingK peaks in the spectrum (61) (69) for the RDE (AE). 2 Refinement: for k = 0 toK − 1 do H k ) by means of (131). end b-Evaluate the kth target amplitude aŝ At the end of the last (namely, the N it th) iteration, an estimate of the complex amplitude of the kth target is evaluated (132)). The ELA can also be employed in the AE; its formulation for the last block can be easily derived from that illustrated above for the RDE by simply replacingĤ If spectral interpolation is used, the term has to be added to C 2D−FFT , in (138); here, I l (I p ) is the number of nodes employed along the first (second) dimension, whereas N l (N p ) is the resulting number of points evaluated by means of interpolation for the first (second) dimension. CSFDEC algorithm -The computational complexity of this method is O(C CSFDEC ), where (see [32, Sec. III-C]) In the last formula, N REF , N it , I M (I N ) are the overall number of re-estimations, the overall number of iterations accomplished in the evaluation of the residuals appearing in (86) and the interpolation order adopted for the first (second) dimension of the considered spectrum in (77), respectively. Alg-P -Unlike the CSFDEC algorithm, this algorithm does not operate in an iterative fashion, performs time domain cancellation and re-computes the spectral residual after each cancellation step in order to get ready for the detection of a new target (if any). For these reasons, its computational complexity is O(C ALG−P ), where In the last formula, C init = M 0 N 0 log 2 (M 0 N 0 ) + M 0 N 0 is the contribution due to step 1) of Alg-P (i.e., to the 2D periodogram maximization), C P = C P M + C P N (with C P M = C P N = 12MN) is the cost originating from the computation of the spectral samples for the first and second frequency of the kth target (see (102) and (103), respectively) and C canc = KMN is the contribution due to cancellation in the time domain. Note that the initialization cost depends on the number of targets K; this is due to the fact that the residual spectrum for the coarse estimation of a new target is evaluated after each time domain cancellation (see (94)). MFA -The initialization phase of this algorithm relies, similarly as both the CSFDEC algorithm and the Alg-P, on the 2D periodogram method; however, since the MFA does not include a cancellation procedure, the search for K local maxima in the periodogram is required in order to acquire K coarse frequency estimates. Moreover, the initialization is followed by a frequency refinement process, which is sequentially repeated for each target. Therefore, the computational complexity of the MFA is O(C MFA ), where In the last formula, C init is equal to C 2D−FFT , in (138), (initialization cost) and QSE -Similarly to the MFA and the ELA, this algorithm evaluates first the coarse estimates of K targets through the 2D periodogram method. The frequency refinement step requires evaluating (120) N it times for each target. Therefore, the computational complexity of the QSE is O(C QSE ), where In the last formula, C init is the initialization cost (which is equal to that of the same step of the MFA and the ELA), N it represents the overall number of iterations carried out to evaluate (120) and C ref = 8MN is the complexity due to the computation of the DFT coefficients required to solve the last referred equation. ELA -The initialization of this algorithm is based on the 2D periodogram method and is followed by the frequency refinement step, which requires solving (131) N it times for each target. Therefore, the computational complexity of the ELA is O(C ELA ), where In the last formula, C init is equal to C 2D−FFT , in (138), (being the cost of 2D periodogram method), N it is the number of iterations carried out to refine the estimates of each target represents the cost of each iteration of the refinement step, being N F D and N F ρ (N F V and N F H ) the sizes of the grid for the refinement of F D and F ρ (F V and F H ), respectively, if the RDE (AE) is considered.

VI. NUMERICAL RESULTS
In our work, seven different embodiments of the DRAEC strategy are compared in terms of computational effort and estimation accuracy achieved in various scenarios. In each embodiment, the same algorithm for the detection and estimation of 2D complex tones is employed in both the RDE and the AE, and the RDE is executed for the second time after that the AE has estimated the DoA (i.e., both the azimuth and the elevation) of all the detected targets in order to generate a finer estimate 19 of their range and velocity. For this reason, in the following, the acronyms FFT0 (FFTi), CSFDEC, Alg-P, MFA, QSE and ELA are adopted to identify the embodiments employing the 2D periodogram method without spectral interpolation (with spectral interpolation), the CSFDEC algorithm, the Alg-P, the MFA, the QSE and the ELA, respectively. It is important to point out that embedding these algorithms in the DRAEC allows us to compare state-of-the-art estimators, in terms of accuracy and complexity, in a 4D radar imaging problem and, in particular, to assess their performance in the estimation of specific parameters of multiple targets. 19. Note that the performance of the second instance of the RDE is affected by the estimation errors introduced by the AE.
In the following, we also assume that: 1) The considered radar system is equipped with a TX HULA (RX VULA) consisting of N T = 8 (N R = 8) elements, whose spacing, as already mentioned in Section II, is d t = λ/2 (d r = λ/2); consequently, the structure of its virtual array is described by Fig. 2. 2) The OFDM modulation employed by the radar system is characterized by the following parameters: a) overall number of subcarriers N = 512; b) overall number of OFDM symbols/frame M = 64; c) subcarrier spacing f = 250 kHz; d) cyclic prefix duration T G = 12.5 μs (consequently, the OFDM symbol duration is T s = 1/ f + T G = 16.5 μs); e) carrier frequency f c = 79 GHz (consequently, the carrier wavelength is λ = c/f c = 3.8 mm); f) cardinality of the PSK constellation N s = 4.
In our simulations, four different scenarios have been considered. The first three scenarios share the following features: 1) They are characterized by a couple of targets (i.e., by K = 2), whose echoes have a unitary amplitude (so that |A 0 | = |A 1 | = 1). The range R 0 , the velocity v 0 , the normalized vertical frequency 20 F V 0 and the normalized horizontal frequency F H 0 of the first target are mutually independent and uniformly distributed random variables (the interval characterizing the uniform distribution of these variables is denoted (X min , X max ) in the following, with X = R, v, F V or F H ). 21 The same parameters for the second target (namely, R 1 ,v 1 , F V 1 and F H 1 ), instead, depend on those of the first one, since they are evaluated as where X = R, v, F V or F H , and X bin and X res represent the normalized target spacing and the resolution, respectively, of the radar system along the X dimension; moreover, R res = c/(2N f ) = 1.1719 m, v res = λ/(2MT s ) = 1.798 m/s, F V res = 1/N R = 0.125 and F H res = 1/N T = 0.125 denote the resolutions in the range, velocity, normalized vertical frequency and normalized horizontal frequency domains, respectively.
2) The signal-to-noise ratio 22 at the RX side varies from −20 to 20 dB. However, the first three scenarios differ for the values selected for the parameters X min , X max and X bin , with X = R, v, F V and F H . In fact, we have that: 1) In the first scenario (denoted S1), ( 20. Note that, given the normalized spatial frequencies F V and F H characterizing a given target, target elevation θ and azimuth φ can be easily computed on the basis of (33) and (34), respectively. 21. In all the scenarios, target parameters have been generated by means of the function rand available in MATLABR2022b. 22. Note that σ 2 W represents the variance of the noise sampleW (p a ,q) m,n by the same parameters in some domain (e.g., by the same DoA).
2) In the first three scenarios, the estimation accuracy achieved by each embodiment of the DRAEC strategy has been assessed by evaluating the root mean square error (RMSE) for the range (X = R), velocity (X = v), azimuth (X = θ) and elevation angle (X = φ) of the considered targets; here, X k [t] denotes the estimate of the parameter X k evaluated for the kth target in the tth Monte Carlo run and N r is the overall number of Monte Carlo runs. Moreover, in applying the last formula, the parameters of the K targets and the estimates generated for them have been ordered on the basis of their range (in particular, according to an ascending order, i.e., from minimum to maximum range).
3) In the fourth scenario, instead, the estimation accuracy achieved by each embodiment of the DRAEC strategy has been assessed by evaluating the normalized RMSE with X = R, v, F V and F H ; here, CRLB X denotes the Cramer-Rao Lower Bound (CRLB) for the estimation of X (the evaluation of the CRLB for the considered scenarios is illustrated in the Appendix-A). Note that the adoption of NRMSE X , in (148), as a performance index allows us to fairly compare the estimation accuracy of each algorithm achieved in the presence of a variable number of targets and its computation is done for each SNR value. 4) The detection thresholds adopted in the inequalities (47) and (52) have not been selected, being K known. In practice, the RDE, in its first instance, searches for K (RDE) = K targets (i.e., K range-Doppler bins). Then, the AE identifies K (AE) ≥ K targets, orders them according to a decreasing perceptual importance and discards the last (K ( In addition, in all the considered scenarios, the following choices have been made for the parameters of the 2D complex tone estimators.
2D periodogram method with interpolation -Orders I l = I p = 7 have been adopted for the spectral interpolation accomplished in both the RDE and the AE, and a grid of size N l × N p = 251 × 251 is selected in the serial refinement of the target estimates.
CSFDEC algorithm -Number of iterations carried out in the evaluation of the residuals N Some numerical results referring to S1 are shown in Fig. 4, where the performance index RMSE X (with X = R, v, θ or φ) characterizing all the considered algorithms is shown for SNR ∈ [ − 20, 20] dB (in these figures and in all the following ones, simulation results are represented by labels, whereas continuous and dotted lines are drawn to ease reading). From these results, it is easily inferred that: 1) The FFT0 is outperformed by all the other methods; note also that the floor observed in the RMSE performance of this embodiment is due to the discretization of the grid (see (76)) employed in the RDE and in the AE (see (61) and (69), respectively).
2) The CSFDEC and the Alg-P achieve very good accuracy (close to the CRLB) thanks to their use of cancellation and refinement procedures.
3) The QSE performs similarly to the Alg-P (FFTi) in range (angle) estimation and similarly to the CSFDEC in velocity estimation.
4) The FFTi takes advantage of peak interpolation, so achieving an estimation accuracy similar to that of the MFA and of the ELA in angular estimation. These considerations, together with those illustrated at point 3), also apply to S2 and S3. 5) The RMSE curves for the FFTi, the MFA and the ELA exhibit a floor at high SNRs. This is due to the fact that the accuracy of the FFTi estimator and the MFA is intrinsically limited by the discretization of their search grid. Further simulation results have evidenced that enlarging the set of trial values improves estimation accuracy; however, this result is achieved at the price of higher computational complexity. As far as the ELA is concerned, its accuracy can also be improved by increasing its number of iterations, but this results in a significant increase in the required computational effort. These considerations apply to all the results shown below for the three algorithms that have been just mentioned.
6) The RMSE R and RMSE v curves of the Alg-P exhibit a floor at high SNRs. This phenomenon can be related to the fact the employed estimation algorithm is biased, since it does not include neither an iterative refinement process nor a leakage compensation procedure for each detected target. 7) The computational efforts required by the FFTi, the CSFDEC, the Alg-P, the MFA, the QSE and the ELA are 1.01, 6.6, 1.1, 10, 1.2 and 49 times higher than that required by the FFT0; these results also hold for S2 and S3.
Further results for S1 are shown in Fig. 5, in which the RMSE R curves of the CSFDEC, the Alg-P and the QSE are shown for the cases in which the RDE is executed only once (dashed lines) and twice (solid lines). From this figure, it is easily inferred that: 1) The improvement in range estimation provided by the second instance of the RDE in order is significant for all the proposed techniques; similar results, not shown here, have been found for velocity estimation.
2) The price to be paid for this improvement is an increase in the overall computational effort (this is quantified by the termN RDE appearing in (137)).
Note also that running a single instance of the RDE corresponds to what is done in all the related technical manuscripts in which only a subset of the target parameters (range and Doppler, in this case) is estimated, whereas the other parameters are kept fixed and/or are not successively compensated for (e.g., see [4] where the angular parameters are not estimated). For this reason, these results evidence the importance of estimating all the target parameters jointly. Some numerical results referring to S2 are shown Fig. 6, where the performance index RMSE X (with X = R, v, θ or φ), characterizing all the considered algorithms, is shown for SNR ∈ [−20, 20] dB. These results lead to the following conclusions: 1) The FFTi, the MFA and the ELA achieve similar accuracy in the estimation of angular parameters and exhibits similar trends (and, in particular, a floor); however, the FFTi is outperformed by the MFA and the ELA in range and Doppler estimation.
2) The RMSE θ and RMSE φ referring to the Alg-P and the QSE are quite flat. This is due to the fact the two targets are located in the same range-Doppler bin. This affects the quality of the signal generated by the AE in compensating for the range and Doppler of each target (see (49)) and passed to the RDE. Moreover, the performance of the QSE is appreciably influenced by the selection of the shifting parameters (i.e., q V and q H for the AE). The values of these parameters have been optimized according to [31,Sec. III,eqs. (40) and (45)].
3) The CSFDEC performs substantially better than all the other embodiments in the estimation of azimuth and elevation, and similarly as the QSE in range and velocity estimation (the accuracy of both embodiments is very close to the CRLB).
Some numerical results obtained for S3 are illustrated in Fig. 7, showing again the dependence of the RMSEs on the SNR, with SNR ∈ [ − 20, 20] dB. These results lead to the following conclusions: 1) The CSFDEC and the QSE achieve the best estimation accuracy (very close to the CRLB) for all the considered parameters, whereas the Alg-P performs similarly in angle estimation only.
2) The accuracy provided by the MFA in azimuth and elevation estimation is slightly better than that characterizing the ELA and the FFTi.
3) The accuracy achieved by the MFA in range and velocity estimation is similar to that provided by the Alg-P and the ELA; moreover, the trend of their RMSE R and RMSE v curves remains flat even at high SNR values.
Our last results refer to S4 and are shown in Figs. 8 and 9. In particular, in Fig. 8 the NRMSE characterizing all the considered embodiments at a given SNR and in the presence of a variable number of targets is shown. These results lead to the following conclusions: 1) The NRMSE X increases with the overall number of targets (i.e., the RMSE X departs from the associated CRLB), with X = R, v, F V or F H . This is due to the fact that increasing K results in a stronger spectral leakage and, consequently, in poorer estimation accuracy of each algorithm.
2) The CSFDEC, the Alg-P and the QSE perform similarly in range and velocity estimation, whereas the CSFDEC performs better in the estimation of azimuth and elevation. This confirms once again that the CSFDEC algorithm better exploits the limited information available in the angular domain 24 and takes advantage of its leakage compensation mechanism. 24. The estimation of angular parameters is based on the N T × N R matrixH k ; see Section III-B. 3) The QSE represents the best option in the case of single target. However, if the overall number of targets increases, it is outperformed by the CSFDEC and the Alg-P; this is due to the fact that the CSFDEC algorithm and the Alg-P make use of a serial cancellation procedure.
4) The ELA performs better than the MFA in the considered scenario.
5) The FFT0 and the FFTi are less accurate than all the other techniques.
In Fig. 9, instead, the computational complexity of all the embodiments is represented for a variable number of targets; both the computational cost, measured in mega FLOPs (MFLOPs), and the computation time 25 (CT) are taken into consideration. From this figure, it is easily inferred that: 1) The trend of most of the CT curves is similar and in agreement with that characterising the corresponding curves of the computational cost; the only exception is represented by the CSFDEC, for which the trend of the computational cost appears to be flatter than that of the CT.
2) The slopes of all the curves (i.e., the relative increase in complexity as K gets larger) are similar.
3) The CT of the Alg-P is very close to that required by the FFT0 and the FFTi; moreover, the last two methods require similar CTs. 4) Even if the complexity of the QSE is similar to that of the Alg-P, the CT of the former embodiment tends to be larger than that of the latter; this is mainly due to the fact their estimators have different initializations. In fact, in the case of the QSE (the Alg-P), the search for K local maxima (for a single maximum) is required.
The results shown in this section evidence that the CSFDEC represents the winning option among the set of considered embodiments, since it achieves the best performance-complexity trade-off.

VII. CONCLUSION
In this manuscript, a novel general strategy to develop suboptimal methods for the detection of multiple targets and the estimation of their parameters in a MIMO OFDM-based JCAS system has been proposed. This strategy is based on the idea of splitting a complicated optimization problem into a couple of simpler (but interacting) sub-problems. Seven different embodiments of it have been described, their complexity has been assessed and their estimation accuracy has been compared in four different scenarios. Our numerical results, based on synthetically generated data referring   to four distinct scenarios, evidence that all the proposed embodiments perform reasonably well, but may require substantially different computational efforts. Moreover, the estimation accuracy of the majority of the algorithms exhibits a floor as the SNR increases. This phenomenon is due to the lack of an iterative procedure for refining the coarse estimates of the detected targets or to the use of a sub-optimal refinement procedure or to the adoption of a (discretized) search grid.
We believe that our work sheds new light on a complicated technical problem, which plays a key role in the development of future JCAS systems; in fact, it provides an in-depth analysis of the accuracy-complexity trade-off characterizing different solutions to it. In a number of applications, achieving good estimation accuracy represents a fundamental requirement; at the same time, real-time operation is also needed, so that substantial attention must be paid to the computational effort required by the adopted estimation algorithms. All in all, we believe the detection and estimation algorithms based on the strategy we propose can represent good candidates for the processing to be accomplished in OFDM-based 4D radars. Our future work includes the application of the proposed strategy to MIMO JCAS systems employing the Orthogonal Time Frequency Space (OTFS) modulation.

A. CRAMER-RAO LOWER BOUND DERIVATION
Cramer-Rao lower bounds for OFDM-based radar systems have been already derived in [4] and [16], but refer to the estimation of Doppler and range only. In this Appendix, the procedure we followed in the evaluation of the CRLBs for Doppler, range, azimuth and elevation is sketched.
First of all, let us consider the signal model (40), that refers to K distinct targets; the parameters of these targets are collected in the vectors 26 γ 0 , γ 1 , . . . , γ K−1 T , (149) and We also define: 1) the trial vectors˜ ,F D ,F ρ ,F V and F H in a similar way as , F D , F ρ , F V and F H , respectively (see (149)-(153)); 2) the vectorsˆ ,F D ,F ρ ,F V and F H , structured like , F D , F ρ , F V and F H , respectively, but collecting the ML estimates of all the considered parameters. The CRLBs we are interested in refer to the ML estimation problem here, is a mean square error (MSE) referring to the (M N N T N R )dimensional column vector and its useful component In the last two formulas, is a (5K)-dimensional column vector collecting the trial values of the target parameters, w is a (M N N T N R )-dimensional noise vector,˜ = (F D ,F ρ ,F V ,F H ), B F ρ k a 0 −F ρ k , a 1 −F ρ k , . . . , a n −F ρ k , . . . , a N−1 −F ρ k T , 26. The kth element of the vector represents the complex amplitude of the kth target and coincides with the parameter A k defined in Section VI. respectively (a z (F X ) is defined by (32), with z = m, n, q or p, and X = D, ρ, V or H if z = m, n, q or p, respectively). If we assume that the elements of the noise vector w are Gaussian, mutually independent and have zero mean and variance σ 2 W , the CRLBs of all the parameters of interest are represented by the diagonal elements of the matrix and ϒ X 0, −j2π, . . . , −j2π(X − 1) T .
for any integer X.
In our work, all the above-mentioned CRLBs have been evaluated numerically on the basis of (164)-(171); in particular, in the computation ofĀ k ,B k ,C k andD k , the following choices have been made for any k: 1) the complex gainγ k (corresponding to A k in Section VI) has been set to unity in all the considered scenarios; 2) the noise variance σ 2 W has been derived from (146), since the SNR and the amplitudes of the target echoes are known in all the considered scenarios. Moreover, as far as the normalized frequenciesF D k , F ρ k ,F V k andF H k are concerned, the following rules have been followed: a)F X 0 has been set to the expected valueF X 0 of the normalized frequency of the first target (with X = D, ρ, V or H). This choice is motivated by the fact that, in all the considered scenarios, the four normalized frequencies characterizing the first target are uniformly distributed random variables; therefore,F X 0 = (F X max − F X min )/2 (the values of F X max and F X min are provided in Section VI).
b) The other (K−1) frequencies {F X k ; k = 1, 2, . . . , K−1} have been set toF X 0 + k F X bin F X res for any k; here, F X bin = X bin represents the normalized bin spacing between adjacent targets adopted in the considered scenario and F X res = 1/Q, where Q = M, N, N R or N T if X = D, ρ, V or H, respectively.