Novel Deterministic Detection and Estimation Algorithms for Colocated Multiple-Input Multiple-Output Radars

In this manuscript, the problem of detecting multiple targets and estimating their spatial coordinates (namely, their range and the direction of arrival of their electromagnetic echoes) in a colocated multiple-input multiple-output radar system operating in a static or slowly changing two-dimensional or three-dimensional propagation scenario is investigated. Various solutions, collectively called range & angle serial cancellation algorithms, are developed for both frequency modulated continuous wave radars and stepped frequency continuous wave radars. Moreover, specific technical problems experienced in their implementation are discussed. Finally, the accuracy achieved by these algorithms in the presence of multiple targets is assessed on the basis of both synthetically generated data and of the measurements acquired through three different multiple-input multiple-output radars and is compared with that provided by other methods based on multidimensional Fourier analysis and multiple signal classification.


I. INTRODUCTION
In the last years, the advances in millimeter-wave semiconductor technology and the development of novel signal processing techniques have lead the way to the exploitation of multiple-input multiple-output (MIMO) radar systems in a number of fields. These systems can be divided in statistical MIMO radars [1], [2] and colocated MIMO radars [3] on the basis of the distance between their transmit array and their receive array; in the first case, transmit and receive antennas are widely separated, whereas, in the second one, they are closely spaced and, in particular, they are usually placed on the same shield. In this paper, we focus on colocated MIMO radars only; such systems play an important role in a number of applications, because of their limited cost, their small size and their ability to detect the presence of multiple targets. The performance achieved by any colocated MIMO radar system depends not only on some important characteristics of its hardware (e.g., the operating frequency, the number of The associate editor coordinating the review of this manuscript and approving it for publication was Jiayi Zhang . transmit and receive antennas, the configuration of the transmit and receive arrays, etc.), but also on the techniques employed in the generation of its radiated waveforms and in the processing of the measurements acquired through its receive array. As far as the last issue is concerned, it is worth stressing that optimal (i.e., maximum likelihood, ML) techniques for the estimation of the overall number of targets and of their spatial coordinates cannot be employed in practice, since they require solving complicated multidimensional optimization problems and, consequently, entail a huge computational effort, even in the presence of a small number of targets [4]. This has motivated the development of various sub-optimal estimation techniques able to achieve good estimation accuracy at a manageable computational cost. A well known sub-optimal technique employed in real world radar systems is the one described in ref. [5] for frequency modulated continuous wave (FMCW) radar systems; this requires: a) the computation of a multidimensional Fast Fourier Transform (FFT) of the matrix collecting the time-domain samples of the signals acquired through the receive array of the employed radar device; b) the search for the peaks of the resulting amplitude spectrum over a range-azimuth-elevation domain or a range-azimuth domain in three-dimensional (3D) and two-dimensional (2D) imaging, respectively. Despite the practical importance of this technique, the computational effort it requires is still significant, since it involves multidimensional spectral analysis of the acquired signals. Moreover, it suffers from the following relevant drawback: it can miss targets whose electromagnetic echoes are weaker than those generated by other spatially close targets; this is due to the fact the spectral contribution due to weak echoes is usually hidden by the leakage originating from stronger echoes. This drawback may substantially affect the overall quality of radar imaging in the presence of extended targets, since such targets can be usually modelled as a cluster of point targets characterized by different radar cross sections [6].
Alternative sub-optimal techniques available in the technical literature are based on the idea of turning a complicated multidimensional optimization problem into a series of simpler and interconnected optimization sub-problems, each of which involves a search for the local maxima of a specific cost function over a limited one-dimensional (1D) or 2D parameter space. Examples of this approach are offered by [7] and [8], and by [9], where range-azimuth estimators for FMCW MIMO radars and for stepped frequency continuous wave (SFCW) MIMO radars, respectively, are derived. More specifically, on the one hand, target delays are first estimated by applying the multiple signal classification (MUSIC) algorithm to a temporal auto-correlation matrix or by identifying the beat frequencies in the downconverted received signal through spectral analysis (in particular, through the FFT algorithm) in [7] and in [8], respectively; then, the acquired information are exploited to estimate the direction of arrival (DOA) of the echoes originating from detected targets. On the other hand, a different approach is proposed in [9], where various iterative deterministic methods applicable to a 2D propagation scenario are derived. These methods have the following relevant features: 1) they process a single snapshot of the received signal (acquired over the whole antenna array); 2) they estimate a new target in each iteration; 3) they do not require prior knowledge of the overall number of targets; 4) they involve 1D or 2D maximizations only; 5) they achieve a good accuracy at a reasonable computational cost; 6) the computational effort they require can be easily controlled by setting a threshold on the maximum number of targets to be detected.
In this manuscript, six new detection and estimation algorithms for 2D and 3D radar imaging are developed. They represent different versions of the same algorithm, called range & angle serial cancellation algorithm (RASCA), and can be interpreted as embodiments of a general approach to target detection and estimation. In addition, they share some important features with the detection and estimation techniques developed in [8] and [9]. In fact, similarly as the techniques illustrated in [9], they are deterministic, process a single snapshot, operate in an iterative fashion and are computationally efficient; the last feature can be related to the fact they require the evaluation of 1D FFTs only and the search for the global maximum of proper cost functions over 1D (frequency, azimuth or elevation) domains. Moreover, similarly as [8], they first extract the range of each detected target from the spectra of the received signals and, then, fuse the information originating from such spectra to extract DOA information. In addition, they exploit the iterative estimation techniques developed in [10] and based on a serial cancellation approach for evaluating the parameters of multiple overlapped sinusoids or multiple overlapped complex exponentials in the presence of additive noise. The devised algorithms are able to accurately detect and estimate multiple close targets, and to solve the problem of merged measurements or unresolved measurements [11]- [14], in the sense that targets generating merged measurements in the range domain are resolved in the estimation of their angular coordinates.
The remaining part of this manuscript is organized as follows. In Section II, the architecture of colocated FMCW/SFCW MIMO radars and the models adopted in our work for representing their received signals are analysed. The general approach to target detection and estimation on which our algorithms are based is illustrated in Section III, whereas the algorithms themselves are described in Sections IV and V, which are devoted to FMCW and SFCW radars, respectively. Various important details about these algorithms are provided in Section VI, whereas some technical limitations encountered in their implementation in real world radar systems are discussed in Section VII. A description of other FFT-based and MUSIC-based detection and estimation algorithms with which our algorithms are compared is provided in Section VIII. The computational cost of all the considered algorithms is illustrated in Section IX, whereas their performance is analysed in Section X, where various numerical results, based on both synthetically generated data and experimental measurements, are illustrated. Finally, some conclusions are offered in Section XI.

II. SIGNAL AND SYSTEM MODELS
This section is devoted to: a) describing the architecture of colocated and bistatic MIMO radar systems operating in time division multiplexing (TDM) [4] and at millimeter waves; b) analysing the model of their received signals in the cases of FMCW and SFCW transmissions. The considered radar systems share the following important characteristics: 1) They are equipped with a transmit (TX) and a receive (RX) antenna array, consisting of N T and N R elements, respectively. 2) These arrays allow to radiate orthogonal waveforms from different antennas and to receive distinct replicas of the electromagnetic echoes generated by multiple targets; moreover, the orthogonality of the transmitted waveforms is achieved by radiating them through distinct TX antennas over disjoint time intervals. 3) Their TX and RX arrays are made of distinct antenna elements, placed at different positions. However,  TX antennas are close to the RX ones and, in particular, are usually placed on the same shield. The architecture of the FMCW and SFCW radar systems which we always refer to in our work are illustrated in Figs. 1a and 1b, respectively. In both cases, the radar transmitter consists of a waveform generator feeding a voltage controlled oscillator (VCO), whose output signal is radiated by the TX array after power amplification. The radiated signal is reflected by multiple targets, whose echoes contribute to the signals acquired through the RX array; each received signal feeds a low noise amplifier (LNA), whose output undergoes downconversion, filtering and analog-to-digital conversion. Finally, the resulting stream of signal samples (i.e., of raw data) is processed for target detection and estimation.
Our derivation of the received signal model for the radar systems shown in Figs. 1a and 1b relies on the following assumptions: a) The radar operates in a static or slowly time varying propagation environment. b) The signal radiated by the radar transmitter is reflected by L static point targets, so that the useful part of the received signal consists of the superposition of L components, each originating from a distinct target. c) All the TX and RX antennas belong to the same planar shield, so that a two-dimensional reference system lying on the physical antenna array can be defined. d) Any couple of physical TX and RX antennas of the considered bistatic radar is replaced by a single virtual antenna of an equivalent monostatic radar.
The abscissa x v and the ordinate y v of the v-th virtual antenna (VA) element associated with the p-th TX antenna and the q-th RX antenna (briefly, the (p, q) VA) are computed as 1 and respectively, with p = 0, 1, . . ., N T − 1, q = 0, 1, . . ., N R − 1 and v = 0, 1, . . ., N VR − 1; here, x p , y p ( x q , y q ) are the coordinates of the p-th TX (q -th RX) antenna and N VR N T · N R represents the overall number of available VAs.

A. MIMO FMCW RADAR SYSTEM
Let us focus now on the FMCW radar shown in Fig. 1a. Its waveform generator produces a periodic sawtooth signal, 1 Note that this is not the only rule adopted in the technical literature to compute the coordinates of the (p, q) VA. For instance, in [15,Par. 4.3.1,, the abscissa (ordinate) of this element is evaluated as 2x v (2y v ), where x v and y v are expressed by (1) and (2), respectively. Keep in mind, however, that if the last rule is adopted, all the following formulas involving such coordinates must be changed accordingly. so that the instantaneous frequency of the chirp frequency modulated signal available at the output of its VCO evolves periodically, as illustrated in Fig. 2a. In this figure, the parameters T , T R and T 0 represent the chirp interval, the reset time and the pulse period (or pulse repetition interval), respectively [4], whereas the parameters f 0 and B are the start frequency and the bandwidth, respectively, of the transmitted signal. Note that, if all the available TX diversity is exploited and a time slot of T 0 s is assigned to each TX antenna, a single transmission frame, over which the transmission from the whole TX array is accomplished, lasts T F N T T 0 s; in this interval, a single snapshot is acquired at the receive side. Let us focus now on a single chirp interval and, in particular, on the time interval (0, T ), and assume that, in that interval, the p-th TX antenna is employed by the considered radar system (with p ∈ {0, 1, . . . , N T − 1}); the signal radiated by that antenna can be expressed as where A RF is its amplitude, and µ B/T is the chirp rate (i.e., the steepness of the generated frequency chirp). It can be proved that, under the assumptions listed above, the n-th received signal sample acquired through the v-th VA element (with v = 0, 1, . . ., N VR − 1) in the considered chirp interval is given by (e.g., see [16,Par. 4.6,eq. (4.27)]) x (v) r,n = L−1 l=0 a l cos 2πn F with n = 0, 1, . . ., N − 1; here, N is the overall number of samples acquired over each chirp period, a l is the amplitude 2 of the l-th component of the useful signal, represents the complex amplitude of the real tone appearing in the right hand side (RHS) of (6), is the normalised version of the frequency characterizing the l-th target detected on the v-th virtual RX antenna, T s (f s ) denotes the sampling period (frequency) of the employed analog-to-digital converters (ADCs), is the delay of the echo generated by the l-th target and observed on the v-th virtual channel, R l , θ l and φ l denote the range of the l-th target, its azimuth and its elevation (all measured with respect to the center of the receive array), respectively, and w r,n is the n-th sample of the AWGN sequence affecting the received signal (this sample is modelled as a real Gaussian random variable having zero mean and variance σ 2 for any v). It is important to point out that: a) the real signal model (6) can be adopted in all the FMCW radar systems acquiring only the in-phase component of the signal captured by each RX antenna; b) some commercial MIMO radar devices provide both the in-phase and quadrature components of the received RF signals (e.g., see [17, Par 2.1 eq. (2.2)]). In the last case, the complex model must be adopted in place of its real counterpart (6) for any n; here, for any v and l, and w (v) c,n is the n-th sample of the AWGN sequence affecting the received signal (this sample is modelled as a complex Gaussian random variable having zero mean and variance σ 2 for any v).

B. MIMO SFCW RADAR SYSTEM
Let us take into consideration now the SFCW radar system shown in Fig. 1b. In this case, the VCO of its transmitter is fed by a staircase generator, so that the instantaneous frequency of the resulting RF signal takes on N f distinct and uniformly spaced values in an interval lasting T s for each TX antenna, as shown in Fig. 2b; the n-th value of the instantaneous frequency is with n = 0, 1, . . ., N f − 1; here, f 0 is the minimum radiated frequency, f is the frequency step size and N f is the overall number of transmitted frequencies.
If the sampling frequency f s = f is assumed for analog-to-digital conversion at the receive side, the measurement acquired through the v-th virtual element at the n -th frequency can be expressed as (e.g., see [9, eq. (3), Sect. II.B]) with v = 0, 1, . . ., N VR − 1, n = 0, 1, . . ., N − 1 and N = N f . Here, and are the complex amplitude and the normalised delay, respectively, characterizing the l-th target and observed on the v-th virtual antenna; moreover, the parameters (a l , τ (v) l ) and the random variable w (v) c,n have exactly the same meaning as the one illustrated for the received signal model described by (12) and the phase ψ (v) l is expressed by (11).

C. TARGET DETECTION AND ESTIMATION
In both the considered FMCW and SFCW radar systems, the useful component of the received signal observed on each virtual channel can be represented as a superposition of L real or complex oscillations; moreover, the value of the parameter L has to be considered unknown. In the following derivations, the real samples {x c,n ; n = 0, 1, . . ., N −1} acquired on the v-th virtual channel are collected in the N -dimensional vector with 3 z = r or c. This vector is processed by the next stages of the radar receiver for target detection and estimation.
As it can be easily inferred from (6)-(9) ((12)-(13) and (15)- (17)), in the considered radar system, the problem of target detection and range estimation on the v-th virtual channel is equivalent to the classic problem of detecting multiple overlapped sinusoids (multiple overlapped complex exponentials) in the presence of AWGN and estimating their frequencies (delays) [18]. In fact, if, in a FMCW radar system, the l -th tone is found at the frequencyf (v) l , the presence of a target at the range (see (9) and (10)) is detected. Similarly, in a SFCW radar system, the normalised delayF (v) l estimated on the v-th virtual channel is associated with a target whose range is (see (10) and (17)) Information about the angular coordinates (namely, the azimuth and the elevation) of this target, instead, can be acquired through the estimation of the set of N VR phases l ; v = 0, 1, . . ., N VR − 1} observed over the available virtual antennas. In fact, since (see (10) and (11)) where λ c f 0 (22) is the wavelength associated with the frequency f 0 , the sequence {ψ (v) l ; v = 0, 1, . . ., N VR − 1} exhibits a periodic behavior characterized by the normalised horizontal spatial frequency if the considered virtual elements form an horizontal uniform linear array (ULA), whose adjacent elements are spaced d H m apart. Dually, if a virtual vertical ULA is assumed, the periodic variations observed in the same sequence of phases are characterized by the normalised vertical spatial frequency where d V denotes the distance between adjacent elements of the virtual array itself. Consequently, angle finding can be easily accomplished by digital beamforming, i.e. by performing a FFT on the estimated phases taken across multiple elements of the virtual array in a single frame interval [19], [20]. Finally, it is important to note that, in the development of detection and estimation algorithms for colocated MIMO radar systems operating at millimeter waves, the following technical issues need to be taken carefully into account: 1) These radar systems often operate at short ranges and in the presence of extended targets. Each radar image is a cloud of point targets whose mutual spacing can be very small [6]. For this reason, the accuracy of these images depends, first of all, on the frequency resolution (delay resolution) achieved by the detection and estimation algorithm employed on each virtual antenna in a FMCW (SFCW) radar system. In fact, this makes the radar receiver able to separate point targets characterized by similar ranges. 2) Distinct radar echoes can be characterized by substantially different signal-to-noise ratios (SNRs), because of relevant differences among the amplitudes of the L overlapped oscillations forming the useful component of the received signal (see (6), (12) and (15)). 3) The number N of samples acquired over each virtual channel usually ranges from few hundreds to few thousands. The last two issues explain why significant attention must be paid to the accuracy achieved by the adopted detection and estimation algorithms at low SNRs and/or for relatively small values of N , since this can appreciably influence the quality of the generated radar image.

III. DESCRIPTION OF THE PROPOSED APPROACH TO THE DETECTION AND ESTIMATION OF MULTIPLE TARGETS
All the algorithms developed in the following section can be considered as specific instances of a general approach to FIGURE 3. Block diagram describing the general approach to target detection and estimation adopted in our work. target detection and estimation; this approach is described by the block diagram shown in Fig. 3. The processing accomplished by the blocks which this diagram consists of can be summarized as follows. Each vector of the set {x (v) z }, collecting N VR vectors (see (18)), undergoes FFT processing, so that, in a FMCW (SFCW) radar system, the analysis of the acquired measurements is moved from the time-domain (frequency-domain) to the frequency-domain (time-domain). The output of the FFT block is processed by the range profile estimator (RPE), that generates the so called target range profile (TRP), i.e. a collection of: a) the ranges at which the relevant echoes are detected; b) the associated energies. Note that the last quantities allow us to rank each range on the basis of its perceptual importance. The output of the FFT processing block and the target range profile are processed by the spatial estimator (SPE). This block detects all the targets associated with each range appearing in the TRP and estimates their angular parameters; moreover, it may generate a finer estimate of their range. The SPE output is represented by the set I t R l ,θ l ,φ l , Ĉ l ; l = 0, 1, . . . ,L − 1 (25) or the set I t R l ,θ l , Ĉ l ; l = 0, 1, . . . ,L − 1 (26) in the case of 3D and 2D imaging, respectively; here, L represents an estimate of the parameter L (i.e., of the overall number of point targets), whereasR l ,θ l ,φ l and |Ĉ l | represent an estimate of the range R l , azimuth θ l , elevation φ l and amplitude |C l |, respectively, of the l-th target (with l = 0, 1, . . .,L − 1). It is important to point out that: 1) If this approach is adopted, range estimation is decoupled from angular estimation, so that a 3D (2D) detection and estimation problem is turned into a) a 1D detection/estimation problem involving the detection of multiple targets and the estimation of their ranges only plus b) a 2D (1D) estimation problem concerning the targets detected at the same range and the estimation of their azimuth and elevation (azimuth only). Consequently, the overall problem of detecting multiple targets and estimating their range and angles is turned into a couple of simpler detection and estimation problems. 2) The SPE exploits the range information generated by the RPE in order to concentrate its computational effort on a set of well defined ranges; this allows to reduce the size of the search space involved in spatial estimation. This explains also why the processing accomplished by the SPE cannot start before that at a least a portion of the range/energy information (i.e., a portion of the TRP) generated by the RPE becomes available. 3) Various techniques can be exploited in the RPE and in the SPE to develop computationally efficient embodiments of the proposed approach.
As far as the last point is concerned, the following techniques can be adopted by the RPE to mitigate its complexity:

a)
Antenna selection -This consists in feeding the RPE with a subset of the outputs of the FFT processing block; such outputs are generated on the basis of N A of the N VR VAs. Note that, on the one hand, a larger N A allows to compute a more accurate TRP; on the other hand, selecting a smaller N A results in a reduction of the overall effort required for the computation of the TRP. b) Antenna-by-antenna processing -The measurements acquired through the selected N A VAs can be efficiently processed by adopting a two-step procedure. In the first step, target range estimation is accomplished on each VA independently of all the other VAs, i.e. the acquired measurements are processed on an antenna-by-antenna basis; this is beneficial when parallel computing hardware is employed in the execution of the first step. In the second step, instead, the target range information extracted from each of the selected N A VAs are fused to generate the TRP. c) Serial target cancellation in the range domain -Target detection and range estimation on each VA represent a multidimensional problem since they aim at detecting multiple targets and estimating their ranges. In our method, this multidimensional problem is turned into a sequence of 2D estimation problems by adopting a serial interference cancellation (SIC) approach (e.g., see [9]). This means that the noisy signal observed on each VA is processed in an iterative fashion. In each iteration, a single (and, in particular, the strongest) target is detected, and its range and complex amplitude are estimated. Then, the contribution of this target to the received signal is estimated and subtracted from the signal itself (i.e., the detected target is treated as a form of interference to be cancelled), so generating a residual signal. The last signal represents the input of the next iteration. This procedure is repeated until the overall energy of the residual drops below a given threshold. Note also that the use of this SIC-based approach allows us to mitigate the impact of the spectral leakage due to strong targets, that can potentially hide weak targets having similar ranges.

d)
Alternating maximization -The estimation of the normalised frequency (or delay) and the complex amplitude of a detected target requires searching for the maximum (or the minimum) of a proper cost function over a 2D domain. In our method, the alternating maximization (AM) technique is exploited to develop iterative algorithms that alternate the estimation of the normalised frequency (or delay) of a given target with that of its complex amplitude; for this reason, a 2D optimization problem is turned into a couple of interacting 1D optimization problems (e.g., see [21, Par. IV-A]).
In the SPE block, instead, the following techniques can be employed to reduce its overall computational complexity:

a)
Alternating maximization -The AM technique is exploited to develop iterative algorithms that alternate the estimation of the elevation of a given target with that of its azimuth. This allows us to decouple the estimation of target elevation from that of its azimuth.

b)
Serial target cancellation in the angular domain -Each of the ranges collected in the TRP is associated with an unknown number of targets; for this reason, the processing accomplished by the SPE aims at resolving all the targets associated with a given range and estimating their angular coordinates.
In the technical literature about radar systems, the detection of an unknown number of targets characterized by the same range (or by ranges whose mutual differences are below the range resolution of the employed radar system) and the estimation of their angular parameters is known to be a difficult multidimensional problem (e.g., see [7, Par. III-C]). In our method, a SIC approach is exploited to turn this multidimensional problem into a sequence of 2D (1D) estimation problems in 3D (2D) imaging (see [9] and references therein). This means that the noisy data referring to a specific range and acquired on all the VAs are iteratively processed to detect a single (and, in particular, the strongest) target, and to estimate its angular coordinates and complex amplitude. Then, the contribution of this target to the outputs of the FFT processing block is estimated and subtracted from them, so generating a set of residual data. This detection/estimation/cancellation procedure is iteratively applied to the residual data until their overall energy drops below a given threshold. Moreover, in a 3D propagation scenario, this procedure is combined with the AM approach described in the previous point; this allows to detect and estimate the angular parameters of a single target (i.e., to solve a 2D optimization problem) by means of an iterative procedure alternating the estimation of its elevation with that of its azimuth (i.e., by means of an algorithm solving a couple of 1D optimization problems). Note also that, once again, the use of a serial cancellation approach allows us to mitigate the impact of the spectral leakage due to strong targets, that can potentially hide weak targets having similar spatial coordinates. c) Parallel processing of the data associated with different ranges -The detection and the estimation of the targets associated with distinct ranges of the TRP can be accomplished in a parallel fashion or in a sequential fashion. The first approach is more efficient than the second one if spatial estimation is executed on parallel computing hardware. In fact, in this case, multiple spatial estimation algorithms can be run in parallel, one for each of the ranges appearing in the TRP. Note, however, that the price to be paid for this is represented by the fact that the target information generated by all the parallel procedures need to be fused when they end. In fact, the analysis of the measurements referring to close ranges appearing in the TRP may lead to detecting the same target more than one time. Based on the general approach outlined above and on the techniques listed for the RPE and the SPE, six specific algorithms, called range & angle serial cancellation algorithms (RASCAs) are developed in the following. Four of these algorithms apply to colocated MIMO FMCW radar systems, whereas the other two to colocated MIMO SFCW radar systems. The algorithms for FMCW radars are called RASCA-FR2 (RASCA-FC2) and RASCA-FR3 (RASCA-FC3), since they generate 2D and 3D radar images, respectively, on the basis of real (complex) measurements. The algorithms for SFCW radars, instead, are dubbed RASCA-S2 and RASCA-S3. In the description of these algorithms we assume, without loosing generality, that the available measurements are acquired through the N VH × N VV virtual uniform rectangular array (URA) represented in Fig. 4 in the case of 3D imaging and through an horizontal ULA (HULA), consisting of N VH virtual antennas, in the case of 2D imaging. In the first case, the horizontal (vertical) spacing between adjacent antennas is denoted d VH (d VV ), whereas, in the second one, is denoted d VH . Moreover, in our considerations, we assume that a reference VA, identified by (p, q) = (p R , q R ) (p = p R ) in the 3D (2D) case, is selected in the virtual array, as exemplified by Fig. 4.
In the following sections, all the RASCAs are described. In Section IV we first focus on the RASCA-FR3 and RASCA-FC3, i.e. on the algorithms to be employed in a MIMO FMCW radar system equipped with a 2D antenna array (in particular, with the URA shown in Fig. 4). Then, we show how to adapt these algorithms to the case in which this radar system is equipped with a 1D antenna array (in particular, with an ULA); this leads to the RASCA-FR2 and the RASCA-FC2. Then, in Section V we show how to adapt RASCA-FC3 and RASCA-FC2 algorithms to be employed in a MIMO SFCW radar system to obtain RASCA-S3 and RASCA-S2, respectively.

IV. RANGE & ANGLE SERIAL CANCELLATION ALGORITHMS FOR A FREQUENCY MODULATED CONTINUOUS WAVE RADAR SYSTEM
In this paragraph, we provide a short description of the architecture of the RASCAs for FMCW radar systems and comment on the method we developed for target detection and cancellation in the angular domain. Then, we illustrate RASCA-FR3 and RASCA-FC3 in detail. Finally, we show how to derive the RASCA-FR2 and RASCA-FC2 from them.

A. ARCHITECTURE OF THE RANGE & ANGLE SERIAL CANCELLATION ALGORITHMS FOR A FREQUENCY MODULATED CONTINUOUS WAVE RADAR SYSTEM
The inner structure of the RASCAs for an FMCW radar system is described by the block diagram shown in Fig. 5, whereas the meaning of the most relevant parameters, sets, vectors and matrices appearing in the description of this algorithm is summarised in Table 1 consisting of 3·N VR N 0 -dimensional vectors. Note that, however, a portion of this set is discarded by the RPE, since this block processes the information originating from N A distinct VAs only. The triads selected by the RPE form the subset of S FFT (27); here,  represents the set of the values of the VA index v identifying the elements of S FFT that belong toS FFT . Each of the triads ofS FFT is processed, independently of all the other ones, by a novel iterative estimation algorithm called single target detection, range estimation and cancellation (STDREC). This algorithm detects the most relevant targets on the selected antenna and estimates their ranges (i.e., the frequencies associated with these ranges; see (9) and (19)) and their complex amplitudes (see (11) and (13)). The name of this algorithm originates from the fact that, in each of its iterations, it detects a single target (namely, the strongest target), estimates its parameters (and, in particular, the frequency characterizing it, i.e. its range) and cancels the target contribution to the received signal; the residual signal resulting from target cancellation represents the input of the next iteration. The output of the STDREC algorithm that processes the raw data originating from the v k -th VA is represented by the set 4 with k = 0, 1, . . ., N A − 1; here, L k is the overall number of targets detected on the considered VA, whereasĈ represent the estimate of the complex amplitude of the i-th target and the index of the frequency bin 5 in which this target has been detected. Finally, the information provided by the N A sets {S v k } are merged to generate the single set where L b is the overall number of targets detected on all the selected VAs,α l is the index of the frequency bin in which the l-th target has been detected and E b,l is the average energy estimated for it. Note that: a) The cardinality L b of the set S RPE represents a preliminary estimate of the overall number of targets; in fact, multiple targets having the same range or ranges whose mutual differences are below the resolution of the employed radar system are detected as a single target and no effort is made at this stage to separate their contributions. b) The energies {E b,l } represent the perceptual importance of the identified frequency bins, in the sense that a larger energy is associated with a more important frequency bin. Both the sets S FFT (27) and S RPE (31) feed the SPE. The aim of this block is to analyse the spectral information associated with the ranges (i.e., with the frequency bins) identified by the RPE in order to: a) estimate the angular coordinates (i.e., azimuth and elevation) of the targets contributing to each frequency bin; b) detect additional targets associated with adjacent frequency bins and potentially hidden by the spectral leakage due to stronger targets; c) estimate the angular coordinates (i.e., azimuth and elevation) of such additional targets and compute a finer estimate of their range.
The first stage of the processing accomplished by the SPE involves the whole set S FFT (27) and is executed on a binby-bin basis, since it aims at: a) detecting all the targets that contribute to the energy of each frequency bin contained in the TRP and b) estimating their angular coordinates. For this reason, this stage consists of L b estimators running in parallel; each estimator focuses on one of the L b frequency bins (i.e., ranges) appearing in the TRP (see Fig. 5). Moreover, each estimator executes a novel iterative estimation algorithm, called single target detection, angular estimation and cancellation (STDAEC). The l-th STDAEC algorithm processes the spectral information available on the whole virtual receive array and referring to theα l -th frequency bin only 4 Note that the complex amplitudeĈ (v k ) i appearing in the following equations is replaced byÂ if the received sequence is complex (see eqs. (7) and (13)). This consideration holds for various equations appearing in the remaining part of this manuscript. 5 Generally speaking, the evaluation of an FFT of order N 0 leads to partitioning the normalised frequency interval [0, 1/2) in N 0 frequency bins.
(with l = 0, 1, . . ., L b −1), detects L[l] targets contributing to it and, for each detected target, computes: a) an estimate of its complex amplitude; b) an estimate of its angular coordinates (i.e., its azimuth and its elevation); c) a refined estimate of its range (do not forget that the preliminary estimate of this range is provided by the bin indexα l ). If D [l] iterations are accomplished by the l-th STDAEC algorithm, D[l] distinct targets are detected in theα l -th frequency bin, provided that none of them is classified as a false (i.e., ghost ) target. In addition, all the estimates generated by this algorithm are collected in the set or in the set Finally, in the second (and last) stage of the SPE, the spatial coordinates of all the detected targets are computed on the basis of the spatial information collected in the L b sets {T l } and an overall image of the propagation scenario is generated in the form of a point cloud.

B. SOME CONSIDERATIONS ON TARGET DETECTION AND CANCELLATION IN THE ANGULAR DOMAIN
It is worth pointing out that the STDAEC algorithm represents the most complicated part of the processing accomplished by all the blocks appearing in Fig. 5. The derivation of this algorithm relies on the fact that: a) each target provides an additive contribution to the spectra evaluated on all the VAs; b) periodic variations are observed in the phase of this contribution if we move horizontally or vertically along the considered virtual array (see Fig. 4). In fact, if we assume that the intensity of the echo received by each VA from the i-th target detected in the l-th frequency bin does not change from antenna to antenna, the complex amplitude C i [p, q, l] observed on the (p, q) VA can be expressed as (see (9) and (10))  here, λ = c/f 0 is the wavelength associated with the start frequency (see (14)), (p R , q R ) is the couple of integers identifying the selected reference VA, θ i [l], φ i [l] and R i [l] are the azimuth, the elevation and the range, respectively, characterizing the considered target, and C i [l] is its complex amplitude observed on the reference antenna. If (34) holds, the rate of the phase variations observed in the complex amplitudes {C i [p, q, l]} for a given l is proportional to (see (23) and (24)) and if we move along an HULA and a vertical ULA (VULA), respectively. In fact, the quantity F H,i [l] (F V,i [l]) represents the normalised horizontal (vertical) spatial frequency characterizing the i-th target detected in the l-th frequency bin; if both these frequencies are known, the elevation and the azimuth of this target can be evaluated as and respectively. Moreover, in the derivation of the STDAEC algorithm, the following two techniques have been exploited: Serial cancellation of targets -This technique is conceptually similar to the cancellation strategy exploited by the STDREC algorithm and allows us to detect multiple targets in the same frequency bin and, in particular, to identify targets having similar spatial coordinates. It is important to keep in mind that the frequencies associated with distinct targets detected in the same frequency bin do not necessarily belongs to that bin; in fact, they can belong to adjacent bins, so that the tails (not the peak) of their spectra are really observed in the considered frequency bin. This problem originates from the fact that, generally speaking, the contribution of a point target to the spectrum computed on each VA is not a line, unless the associated normalised frequency is exactly a multiple of the fundamental frequency consequently, such a contribution is spread over multiple adjacent frequency bins (i.e., spectral leakage is observed) Spatial folding -As already stated above, the frequency associated with a target detected in the l-th frequency bin does not necessarily fall exactly in that bin. The technique dubbed spatial folding has been devised to: a) evaluate a more accurate estimate of the frequency associated with a target detected in a given bin; b) discriminate real targets from ghost targets. Spatial folding is based on the following idea. Once the horizontal and the vertical spatial frequencies associated with a target detected in a given frequency bin have been estimated (see (35) and (36)), the spectra computed on multiple VAs can be combined in a constructive fashion by 1) taking a reference VA (identified by (p, q) = (p R , q R ); see Fig. 4), and compensating for the phase differences, estimated for that target, between the reference VA and the other VAs of the whole array, or 2) taking a reference ULA and compensating for the phase differences, estimated for that target, between the reference ULA and other ULAs parallel to it.
In case 1), folding generates a single spectrum, dubbed folded spectrum , and has the beneficial effects of a) averaging out the effects of the noise that affects the VAs and b) combining, in a constructive fashion, the contributions of all the targets different from the one which the employed spatial frequencies refer to. For this reason, in analysing the amplitude of the folded spectrum, a well defined peak in its amplitude is expected in the l-th frequency bin or in a bin close to it. When this peak is detected, a refined estimate of the frequency (and, consequently, of the range) and the complex amplitude characterizing the target for which folding has been accomplished can be computed by identifying its position. On the contrary, if no peak is found, the detected target is actually a ghost target. In case 2), folding generates as many folded spectra as the number of antennas of the reference ULA and offers the same advantages as case 1).
In the remaining part of this manuscript, when folding is employed, the following terminology is adopted: Vertical folding -This refers to the case in which folding involves a reference HULA on which other HULAs are folded.
Overall folding -This refers to the case in which folding involves all the spectra, i.e. the overall virtual URA; a single folded spectrum is computed in this case.
Note that, in any case, folding may involve the whole virtual receive array or a portion of it. The exploitation of a subset of the available VAs is motivated by the fact that, in practice, in computing a folded spectrum that refers to the l-th frequency bin, the estimatesF H,i [l] andF V,i [l] of the frequencies F H,i [l] and F V,i [l], respectively, are employed, so that the quality of the phase compensation factors computed for the antennas that are farther from the reference antenna or the reference HULA may be affected by significant estimation errors.
All the mathematical details about vertical and overall folding can be found in the next paragraph.

C. DETAILED DESCRIPTION OF THE RANGE & ANGLE SERIAL CANCELLATION ALGORITHMS FOR A FREQUENCY MODULATED CONTINUOUS WAVE RADAR SYSTEM ENDOWED WITH A TWO-DIMENSIONAL ANTENNA ARRAY
In the following, the RASCA-FR3 is described first; then, the (minor) modifications required to obtain the RASCA-FC3 from it are illustrated. The RASCA-FR3 processing is divided in three tasks, each associated with one of the blocks appearing in Fig. 5 (the i-th task is denoted Ti); a description of each task is provided below. Various details about the techniques employed in these tasks, omitted here to ease the understanding of the overall flow of the algorithm, are provided in Section VI.

1) T1-FFT PROCESSING
The processing accomplished within this task can be summarized as follows. Given the vector x with n = 0, 1, . . ., N − 1 and m = 1, 2. Then, the vectors x respectively; here, M is a positive integer (dubbed oversampling factor), 0 D is a D-dimensional (column) null vector and Finally, the N 0 -dimensional vectors with m = 0, 1, 2, are computed for any v (i.e., for any p and q) by executing a N 0 order FFT for each of them; here, DFT N 0 [x] denotes, up to a scale factor, the N 0 order discrete Fourier transform (DFT) of the N 0 -dimensional vector x.

2) T2-RPE
The processing accomplished within this task consists of the three consecutive steps listed below (the i-th step is denoted T2-Si in the following); each step is associated with one of the blocks included in the RPE, as shown in Fig. 5. T2-S1) VA selection -In this step, the setS FFT (28) is built. This requires generating the set S A (29), i.e. a set of N A integers that identifies the selected VAs. In our computer simulations, the elements of S A have been generated by randomly extracting N A distinct integers from the set {0, 1, . . . , T2-S2) Target detection and range estimation -The processing carried out within this step is executed by the STDREC algorithm; this operates on an antenna-by-antenna basis. The STDREC processing for the v k -th VA (with k = 0, 1, . . ., N A − 1) can be summarized as follows. A simple initialization is accomplished first by setting with m = 0, 1, 2, and the iteration index i to 0. Then, the STDREC iterations are started; in the i-th iteration, the three steps described below are accomplished to detect a new target and cancel its contribution to the triad (X STDREC-S1) Detection of a new target and estimation of its parameters -The triad (X is processed to detect a new (i.e., the i-th) target, and to estimate the normalised frequency F is not a multiple of the fundamental frequency F FFT (39), that characterizes the FFT processing executed in T1; for this reason, it can be expressed as where F is a real parameter called residual. This step consists in executing an algorithm, dubbed single frequency estimator 6 (SFE) and whose detailed description is provided in Paragraph VI-A. In short, the SFE computes the estimatesĈ , respectively, on the basis of the triad (X represents the index of the frequency bin in which the i-th target is detected on the v k -th antenna. Note that the parameterF , even if useless in the construction of the set S v k (30), is exploited in the next step.

STDREC-S2) Cancellation of the new target
, given by the i-th (i.e., by the last) target detected on the v k -th VA, to the triad (X is computed on the basis of (130)-(132) (see Paragraph VI-D) and cancelled from the triad itself. Cancellation consists in the computation of the new residual triad with m = 0, 1, 2. 6 Note that our general description of the SFE includes the computation of three DFTs, that, in this case, are already evaluated in T1. 53) is computed and compared with the positive threshold T STDREC (which may depend on range, i.e. on the detected frequency). If this energy is below the threshold, the STDREC algorithm stops and L k = i relevant targets are detected on the v k -th VA; otherwise, the recursion index i is increased by one and a new recursion is started by going back to STDREC-S1. T2-S3) Fusion of range information -This step aims at merging the information provided by the N A sets {S v k } (30) evaluated in the previous step. Its output is represented by the set S RPE (31), whose elements (i.e., the L b couples

STDREC-S3) Computation of the residual energy in the frequency domain -The energy
identifying all the bins in which at least one target has been detected on the v k -th VA (with k = 0, 1, . . ., N A − 1), the set is generated by putting together all the distinct integers that appear at least once in the N A sets {A where represents the overall number of antennas that contribute to this energy (here, The processing accomplished within this task consists of the two consecutive steps listed below (the i-th step is denoted T3-Si in the following); each step is associated with one of the blocks contained in the SPE represented in Fig. 5. T3-S1) Bin analysis -Within this step, L b STDAEC algorithms are run in parallel, one for each of the L b ranges (i.e., frequency bins) appearing in the TRP. A schematic description of l -th STDAEC algorithm is provided below (with l = 0, 1, . . ., L b − 1). This algorithm consists of three steps (its r-th step is denoted STDAEC-Sr in the following) and is initialised by 1) Setting its iteration index i to 0.
2) Setting where is a N VH × N VV matrix collecting the spectral information available on the whole virtual receive array and referring to theα l -th frequency bin only. Then, the STDAEC algorithm starts executing its iterations. Within its i -th iteration, it accomplishes the three steps described below.
STDAEC-S1) Detection of a new target and estimation of its angular parameters -In this step, the N VH × N VV matrix is processed to detect the strongest target contributing to it, and to compute the , respectively (note that this target represents the i-th one detected in the considered frequency bin, since (i − 1) targets have been detected in the previous recursions). This result is achieved by executing a novel iterative detection and estimation algorithm called single target detection and angular estimation (STDAE), whose description is provided after illustrating the overall structure of the RASCA-FR3 to ease reading.

STDAEC-S3) Residual energy test -The energy
of the residual spectrum vector X (i+1) [l] (62) is compared with the positive threshold T STDAEC (which may depend on angular coordinates). If this energy is below the threshold, the STDAEC algorithm stops; otherwise, the recursion index i is increased by one and a new iteration is started by going back to STDAEC-S1. If D [l] iterations are accomplished by the STDAEC algorithm operating on theα l -th frequency bin, no more than D[l] distinct targets are identified in that bin (D[l] targets are found if none of them is deemed to be a ghost target). All the targets information acquired from theα l -th frequency bin are collected in the set T l (32). T3-S2) Evaluation of the target spatial coordinates and generation of the overall image -In this step, the estimates of the range, of the elevation and of the azimuth of the i-th target detected in theα l -th frequency bin are computed as (see (19), (37) and (38) (8)). Finally, these information are fused to generate the overall set I t (25), describing the generated radar image; in general, this image is a cloud ofL points. The set I t (25) results from the union of all the sets {I This concludes our description of the RASCA-FR3. Let us focus now on the most complicated part of the STDAEC algorithm, i.e. on the STDAE algorithm. This algorithm makes use of the so called spatial folding (see the previous paragraph). The exploitation of this procedure in the STDAE algorithm requires: 1) Selecting a reference VULA, that consists of N VULA adjacent and vertically aligned VAs (with N VULA ≤ N VV ), within the virtual array; in the following, we assume, without any loss of generality, that the reference VULA includes the reference antenna and, consequently, is identified by p = p R and q = q I , q I + 1, . . ., q F (with q I ≤ q R ≤ q F ), so that N VULA = q F − q I + 1 (see Fig. 6a). 2) Selecting a reference HULA, that consists of N HULA adjacent and horizontally aligned VAs; in the following, we assume, without any loss of generality, that the reference HULA is the horizontal ULA containing the reference antenna and, consequently, is identified by p = p I , p I + 1, . . ., p F (with p I ≤ p R ≤ p F ) and q = q R , so that N HULA = p F − p I + 1 (see Fig. 6a). 3) Selecting a set of HULAs, different from the reference HULA and having the same size of it (i.e., the same number of VAs); in the following, we assume, without any loss of generality, that these HULAs, called vertically folded HULAs, correspond to q = q Fig. 6b; note that the overall number of involved HULAs is N The STDAE algorithm consists in the four steps described below (its r-th step is denoted STDAE-Sr in the following). STDAE-S1) FFT processing on the reference VULA and vertical frequency estimation -The portion of the initial spectral information referring to the reference VULA is extracted from the matrix X (i) [l] and stored in the N VULA -dimensional vector  that is processed by the complex single frequency estimator 7 (CSFE). This algorithm detects the i-th (strongest target) appearing inα l -th frequency bin and computes the esti- (36)), respectively. Note that the quantityĈ V,i [l] is not exploited in the following since, it represents a preliminary estimate of C i [l].
It is worth pointing out that the execution of the CSFE entails: a) The evaluation of the N VULA -dimensional vector with k = 1 and 2; here, b) The computation of anN 0 order FFT of the vector S andM represents the adopted oversampling factor. This produces the vector with k = 0, 1 and 2. Note that the m-th element of the vector s can be expressed as with m = 0, 1, . . .,N 0 − 1.

STDAE-S2)
Vertical folding -The estimateF V,i [l] of the normalised vertical frequency F V,i [l] (36) is employed to compensate for the phase difference between each of the HULAs selected for vertical folding and the reference HULA (i.e., for the phase differences along the vertical direction), so that the spectral information associated with all these HULAs can be combined (i.e., summed) in a constructive fashion. To this aim, the phase rotation factor is computed for the q-th HULA, with q = q F . Then, vertical folding is accomplished by computing the N HULA -dimensional vector that collects the values taken on by the N HULA vertically folded spectra referring to theα l -th frequency bin; here, is a N HULA -dimensional row vector extracted from the q-th row of the matrix X (i) [l] (61). STDAE-S3) FFT processing and horizontal frequency estimation -The processing accomplished in this step is very similar to that carried out in STDAE-S1. In fact, the only difference is represented by the fact that the STDAE-S4) Overall folding and frequency/amplitude estimation -In this step, the angular information i.e., the fre-quenciesF V,i [l] andF H,i [l] computed in STDAE-S2 and STDAE-S3, respectively, are exploited to accomplish overall folding 8 ; this step involves the whole spectrum computed on the selected VAs. If the whole receive antenna array is exploited, overall folding consists in computing the N 0 -dimensional vector is a phase rotation factor, R and 78), the sequence of the absolute values of its elements is analysed to verify the presence of a peak in thê α l -th frequency bin or in a bin close to it. To this aim, after evaluatingα OF arg max the quantity dα[l] |α OF −α l | is compared with the positive threshold T OF . If dα[l] exceeds T OF , the presence of a ghost target is detected; otherwise, the N 0 -dimensional vector is computed for m = 1 and 2, and the CSFE algorithm 9 is run to estimate, on the basis of the triad (X 0, , respectively; these parameters characterize the i-th target detected in theα l -th frequency bin. Note that the integer part (see (51)) does not necessarily coincide withα l but, if it differs, it is certainly close toα l . Ifα i [l] is different fromα l and appears in one of the couples forming the set S RPE (31), it is discarded, because the corresponding frequency bin is already being analysed by one of the other STDAEC algorithms. Otherwise, the new couple where , is added to the set S RPE and the number of its elements (i.e., L b ) is increased by one. This means that an additional STDAEC algorithm is run on the (new)α i [l]-th bin. This concludes our description of the STDAE algorithm and, consequently, of the RASCA-FR3, whose overall structure is summarised in Algorithm 1.

4) ADDITIONAL COMMENTS
The structure of the RASCA-FR3 (RASCA-FC3) deserves a number of comments, that are listed below for the different tasks and the steps they consist of. T1 -In this task, each of the vectors {X  47) and (18)). The vectors X (v) 1 and X (v) 2 , instead, consist of, up to a scale factor, N 0 equally spaced samples of the first and the second derivatives, respectively, of the same spectrum.
T2-S1 -The exploitation of a subset of the available antennas is motivated by the need of reducing the computational effort required by T2 as much as possible. The adoption of a deterministic method for the selection of N A antennas (with N A < N VR ) is not recommended. In fact, when multiple consecutive snapshots are processed to generate independent images, randomly changing the subset of N A antennas from snapshot to snapshot allows the considered radar system to benefit from antenna diversity.
T2-S2 -The STDREC algorithm deserves the following comments: a) The availability of accurate estimates of the normalised frequency F (7) and (13)) plays an important role in this step, since these parameters are exploited in the serial cancellation procedure based on (53). In particular, ignoring the frequency residual δ 2 } according to (47). 2 T2 -RPE: S1) Extract N A VAs from all the available VAs; then, build the setS FFT (28). (49)); then, set the iteration index i = 0 and compute the initial energy E A threshold on the maximum computational effort required by the STDREC algorithm can be set by requiring that the recursion index i never exceeds a fixed threshold; this is equivalent to limit the overall number of targets that can be detected on each VA. c) The STDREC algorithm generates N A different data sets; the k -th data set consists of the triads {(α , characterizing the L k targets detected on the v k -th antenna (with k = 0, 1, . . ., N A − 1). Note that the overall number of targets may change from antenna to antenna, especially in the presence of extended targets; this is due to the fact that the signals acquired on different VAs can exhibit significant differences in their spectral content. d) The following important interpretation of the processing accomplished by the STDREC algorithm on the v k -th VA can be given. The vector X can be seen as a collection of noisy spectral information referring to N 0 distinct frequency bins (i.e., to N 0 distinct range bins) and is usually dense in the presence of multiple extended targets, as illustrated in Fig. 7-a) (where the absolute value of its elements is represented). The STDREC allows to extract a discrete frequency (i.e., range) profile from the vector X (v k ) 0 , as illustrated in Fig. 7-b). In various real world scenarios, this profile turns out to be sparse, even in the presence of a dense vector X (v k ) 0 ; this is beneficial, since allows to concentrate the RPE computational effort on a set of specific ranges (i.e., frequency bins). The range profile characterizing the v k -th VA is described by the set of L k identifies the frequency bin associated with the i-th target detected on the considered VA, whereas the absolute value ofĈ i ) represents an estimate of the strength of the echo associated with it. e) The STDREC algorithm can be used for detecting multiple targets and accurately estimating their range in a monostatic radar. f) The STDREC algorithm can be easily extended in a way that multiple targets are detected and estimated in parallel in each of its iterations. If we focus on its i-iteration and the v k -th VA, this result is achieved by running multiple (say, m i ) instances of the SFE (CSFE) algorithm in parallel. Each of these instances is initialised with the frequency corresponding to the absolute maximum or a relative maximum detected in the sequence of the absolute values of the elements of the vector X (47)). In this case, a constraint is set on the minimum spacing between the m (v k ) i detected frequencies in order to minimize the interference between the instances running in parallel. Moreover, after identifying the absolute maximum in the above mentioned sequence, a threshold, proportional to such a maximum, is set on the minimum value of the acceptable relative maximum/maxima, so that unrelevant frequencies are discarded. It is also worth stressing that, if a cluster of m ) appearing in the RHS of (53) consists of the sum of m (v k ) i terms, each associated with one of these frequencies. g) The STDREC algorithm employed in the RASCA-FR3 (RASCA-FC3) represents an instance of the single frequency estimation and cancellation (complex single frequency estimation and cancellation) algorithm derived in [10] for the estimation of multiple overlapped real (complex) tones. For this reason, in the case of complex received signals, it can be replaced by one of the multiple tone estimators available in the technical literature, like the CFH algorithm [22], the algorithm developed by Ye and Aboutanios [23], [24] and the algorithm derived by Serbes [25] (the last two algorithms are denoted Alg-YA and Alg-S, respectively, in the following). In fact, all these algorithms are recursive and rely on a serial cancellation procedure since, within each recursion, they detect a single tone, estimate its parameters and subtract its contribution from the residual signal emerging from the previous iteration. h) The estimates generated by the STDREC algorithm are potentially biased if the parameters of the SFE (CSFE) executed in its first step are not properly selected (see [10]). In principle, this bias can be arbitrarily reduced by increasing the overall number of iterations and/or re-estimations accomplished by the SFE (CSFE). However, we found out that, in the case of complex received signal, a computationally efficient alternative to this approach is represented by running an additional step (i.e., STDREC-S4) after that the first three steps of the STDREC algorithm has been carried out. In this final step, the Alg-YA is run after initializing it with the estimates of the normalised frequencies and the associated complex amplitudes generated by the STDREC. The hybrid technique that results from interconnecting the STDREC algorithm with the above mentioned algorithm is dubbed hybrid STDREC (HSTDREC) in the following; note that this algorithm represents an instance of the hybrid CSFE proposed in [10].
T3-S1 -This step is the most complicated of the whole algorithm and deserves the following comments: a) In principle, the horizontal and vertical spatial frequencies (see (35) and (36)) of multiple targets contributing to theα l -th frequency bin can be detected by first computing a 2D DFT of the matrix X [l] (60) and, then, by looking for local maxima over the absolute values of the elements of the resulting 2D matrix; note that the matrix X [l] can be also zero-padded before computing its 2D FFT to improve the resulting spectral resolution. This procedure may require a significant compu-VOLUME 10, 2022 tational effort and its accuracy is affected by the spectral leakage due to any potential strong target.
In the STDAEC algorithm, instead, 2D processing is avoided by alternating vertical and horizontal 1D FFTs. Consequently, relevant spatial frequencies are estimated by searching for the peaks of 1D amplitude spectra (i.e., in the absolute values of the elements of the vectors S [l]); in other words, an AM approach is adopted. Note that this approach allow us to mitigate the overall computational complexity and to detect weak targets hidden by close strong targets through successive cancellations. b) In STDAE-S1, each of the three vectors {s c) The processing accomplished in STDAE-S3 is very similar to that carried out in STDAE-S1.
In fact, the only difference is represented by the fact that the N VULA -dimensional vector S Similarly as the STDREC algorithm, the STDAEC algorithm can also be considered as an instance of the CSFEC algorithm mentioned at point g) of T2-S2. Therefore, in principle, it can be replaced by the CFH algorithm [22], the Alg-YA [23], [24] or the Alg-S [25]. Moreover, a further (and final) step, based the Alg-YA can be added to the STDAEC algorithm to mitigate its estimation bias. e) As already suggested for the STDREC algorithm, the STDAEC algorithm can be employed to detect and estimate multiple angles in parallel; this requires running multiple instances of the CSFE algorithm in parallel.
Our final comments concern the use of RASCA-FR3 and RASCA-FC3 in FMCW radar systems whose virtual antenna array is not an URA; for instance, in our experimental work (see Section X), a colocated MIMO FMCW radar equipped with the virtual receive array shown in Fig. 8 has been employed. Note that the first two processing tasks of the RASCAs are carried out on an antenna-by-antenna basis; therefore, they are not influenced by the shape of the con- sidered virtual array. However, this shape influences the way spatial folding is accomplished in T3. More specifically, as far as the last point is concerned, the following considerations can be made: 1) The array structure represented in Fig. 8 can be treated as an URA if its gaps are zero-padded.
2) The reference VULA should be selected in a way to maximize the number of non-zero vertically aligned VAs and, consequently, the number of VAs contributing to the estimation of the elevation angle, as illustrated in Fig. 8.
3) The reference HULA should be selected in the middle of the antenna array; this mitigates the effects of the errors affecting the estimate of normalised vertical frequency in the vertical folding procedure (do not forget that such errors may have a significant impact on the contributions of the HULAs that are farther from the reference HULA; see (74)).
4) The vertical folding accomplished by the STDAE algorithm involves VULAs of different sizes. More specifically, in the i-th iteration of the STDAEC algorithm, vertical folding is accomplished by computing the N HULA -dimensional vector (see (61) and STDAE-S2) The algorithms described in the previous paragraph can be easily adapted to the case in which the considered colocated MIMO radar system is equipped with a single ULA and, consequently, can be exploited for 2D imaging only; this leads to RASCA-FR2 and RASCA-FC2. The changes made in RASCA-FR3 and RASCA-FC3 to obtain RASCA-FR2 and RASCA-FC2, respectively, concern only the SPE and can be summarized as follows: 1) The first three steps of the STDAE in T3-S1 are not performed; therefore, the fourth step of that algorithm is the first one to be executed. Moreover, the matrix X (i) [l] (61) is replaced by the N VH -dimensional vector collecting the spectral information available on the whole virtual receive array and referring to theα l -th frequency bin only.
2) The spatial frequencyF V,i [l] is unavailable and, therefore, it is not included in the set T l (32); note that the elevation angleφ i [l] (65) is not estimated in this case.

V. RANGE & ANGLE SERIAL CANCELLATION ALGORITHMS FOR A COLOCATED MIMO SFCW RADAR SYSTEM
In this paragraph we analyse the main changes that need to be introduced in the RASCA-FC2 and RASCA-FC3 described in the previous paragraph to adapt them to a MIMO SFCW radar system, i.e. to develop RASCA-S2 and RASCA-S3. Let us focus on RASCA-S3 first. As far as T1 is concerned, the only required change is represented by the replacement of: a) the sequence {x The changes to be made in T2 concern its second step (i.e., the STDREC algorithm), since the noisy measurements processed in this step are always complex (see (15) and (17)); the new structure of this step is summarized below.
T2-S2) Target detection and range estimation -In this step, that aims at detecting the most relevant targets on each of the N A antennas, a minor change is required in the cancellation procedure with respect to its counterpart employed in the MIMO FMCW radar system. This is due to the fact that, as already stated above, the noisy measurements processed in a MIMO SFCW radar system are always represented by complex (frequency-domain) sequences.
The initialization of the STDREC algorithm remains unchanged: the triad (X ) is defined according to (49) (with k = 0, 1, . . ., N A − 1). In its i-th iteration, this algorithm accomplishes the three steps described below (these are denoted STDREC-Sp in the following, with p = 1, 2 and 3).

STDREC-S1) Detection of a new target and estimation of its parameters -In this step, the normalised delay F
(v k ) i and the complex amplitude A (v k ) i associated with the i-th target are estimated. Note that, generally speaking, the normalised delay F is not a multiple of the fundamental delay that characterizes the output of the IFFT processing executed in task T1; for this reason, it can be expressed as where F c,i is a given coarse estimate of the normalised delay F is a real parameter called residual. The processing accomplished in this step is executed by an algorithm dubbed complex single delay estimator (CSDE) in the following (see Paragraph VI-C); this computes the estimatesÂ , respectively, on the basis of the triad (X

STDREC-S2) Cancellation of the new target -The contributions (C
given by the i -th (i.e., by the last) target detected on the v k -th VA to the triad (X , are computed on the basis of (146)-(148) (see Paragraph VI-E). Therefore, the cancellation procedure consists again in the computation of the new residual triad on the basis of (53).

STDREC-S3) Computation of the residual energy in the time domain -The energy of the residual time-domain vector X
is computed (see (54)) and compared with the positive threshold T STDREC . If this energy is below the threshold, the STDREC algorithm stops and L k = i relevant targets are detected on the v k -th VA; otherwise, the recursion index i is increased by one and a new recursion is started (i.e., we go back to STDREC-S1).
T2-S3) Range information fusion -In this case, the set A b (see (56)) collects L b relevant time bins {α l ; l = 0, 1, . . . , L b − 1}, whereas the set {E b,l } contains the energy associated with each of them (see (57)). The average energy E b,l associated with theα l -th delay bin is evaluated as with l = 0, 1, . . . , L b − 1; here, N b,l is defined by (58).
The final output of the information fusion is represented by the set S RPE (31). Note that the cardinality L b of the last set represents a preliminary estimate of the overall number of targets. The processing accomplished in T3 (i.e., by the SPE) has the same structure and interpretation as that illustrated for the FMCW radar system since it aims at estimating the angular VOLUME 10, 2022 parameters of the targets on the basis of the discrete range profile computed in T2. However, frequency bins are now replaced by delay bins. Note also that, in step T3-S1, the matrix X[l] is still defined on the basis of (60), but the complex amplitude A i [p, q, l] appearing in (34) is replaced by

STDAEC-S3)
Residual energy test -In this step, the energy E (i+1) [l] of the residual X (i+1) [l] evaluated in the previous step is computed on the basis of (63)) and is compared with a positive threshold; if this energy is below the threshold, the STDAEC algorithm stops, otherwise the recursion index i is increased by one and a new iteration is started by going back to STDAEC-S1.
All the target information acquired from theα l -th frequency bin are collected in the set is evaluated by the CSDE algorithm (see Paragraph VI-C), for k = 0, 1 and 2; note that with m = 0, 1, . . .,N 0 − 1, and that the quantity S STDAE-S2) Vertical folding -Vertical folding is employed to compensate for the phase differences between the considered HULAs. However, the phase rotation factor appearing in (76) is computed as where andF V,i [l] is the estimate of the normalised vertical frequency F V,i [l] evaluated in STDAE-S1.   (85)) is added to the set S RPE and the number of its elements (i.e., L b ) is increased by one. The RASCA-S3 can be easily adapted to the case in which the radar system is equipped with a single ULA. The changes to be made are exactly the same as those illustrated in the previous paragraph for deriving RASCA-FR2 (RASCA-FC2) from RASCA-FR3 (RASCA-FC3).

VI. DESCRIPTION OF VARIOUS ALGORITHMS EMPLOYED IN THE PROPOSED EMBODIMENTS
In this section, various mathematical details about the techniques employed in the RASCAs are provided.
Since the processing accomplished in T1 of the RASCAs has been fully analysed in the previous section, in this paragraph we provide a detailed description of: a) the SFE (see T2-S2); b) the CSFE (see T2-S2 and T3-S1); c) the CSDE (see T2-S2 and T3-S1); d) the target cancellation procedures employed in T2-S3 and T3-S1.

A. SINGLE FREQUENCY ESTIMATOR
In this Paragraph, the SFE derived in [10] is summarized. This algorithm processes the samples of the real sequence {x r,n ; n = 0, 1, . . ., N − 1}, whose n-th element is with n = 0, 1, . . ., N − 1, and generates an estimate of the normalised frequency F and of the complex amplitude of the real tone appearing in the RHS of (102); here, N is the overall number of elements of the sequence {x r,n }, a and ψ are the tone amplitude and phase, respectively, and {w r,n ; n = 0, 1, . . ., N − 1} is a real AWGN sequence. This algorithm is initialised by 1) Evaluating: a) the vector 0,ZP , (104) where the DFT order N 0 is defined by (46), M is the oversampling factor and x 0 x r,0 , x r,1 , . . . , x r,N −1 where the integerα is computed aŝ c) the quantityρ (0) d) the initial estimateĈ (0) of C aŝ where and respectively; f) the initial estimateˆ (0) of aŝ where g) the first fine estimateF (0) of F aŝ 2) Setting its iteration index i to 1. Then, an iterative procedure is started. The i-th iteration is fed by the estimatesF (i−1) andĈ (i−1) of F and C, respectively, and produces the new estimatesF (i) andĈ (i) of the same quantities (with i = 1, 2, . . ., N SFE , where N SFE represents the overall number of iterations); the procedure employed for the evaluation ofF (i) andĈ (i) consists of the two steps described below (the p-th step is denoted SFE-Sp).
SFE-S1) -The new estimateˆ (i) of is computed as 10 (see (118)-(119)) in the evaluation of the coefficients {b(ρ), c(ρ)} appearing in the RHS of (119),Ĉ =Ĉ (i−1) and 10 The quantities {X k,ρ ; k = 1, 2} required in the computation of the coefficients b(ρ) and c(ρ) can be also evaluated by means of the interpolationbased method illustrated in [10, Sect. III, p. 12]. In our work, barycentric interpolation has been always used [26]; in the following, the parameter I represents the interpolation order. These considerations hold also for the CSFE and the CSDE described below. VOLUME 10, 2022 are assumed. Then, is evaluated.
At the end of the last (i.e., of the N SFE -th) iteration, the fine estimatesF =F (N SFE ) andĈ =Ĉ (N SFE ) of F and C, respectively, become available.

B. COMPLEX SINGLE FREQUENCY ESTIMATOR
All the results illustrated in the previous paragraph refer to the real sequence {x r,n }, whose n-th element is expressed by (102). However, a similar estimation method (namely, the CSFE) has been developed for the complex counterpart, i.e. for a complex sequence {x c,n ; n = 0, 1, . . ., N − 1}, whose n-th element is respectively.

C. COMPLEX SINGLE DELAY ESTIMATOR
The SFE (CSFE) described in Paragraph VI-A (VI-B) is fed by the real (complex) sequence {x r,n } ({x c,n }), whose n-th element is expressed by (102) ((124)); as shown in Section II, these signal models are suitable to represent the sampled downconverted signal available at the receive side of a colocated FMCW MIMO radar in the presence of a single target (i.e., when L = 1). However, a similar estimation method, called CSDE, can be developed for a colocated SFCW MIMO radar system, since the complex (frequency-domain) signal model (15) characterizing this system is very similar to the complex (time-domain) signal model (12) illustrated for a FMCW radar system, the only differences being represented by: a) the physical meaning of the parameter F (v) l , since this represents a normalised delay in (15) and a normalised frequency in (12)); b) the sign of the argument of all the complex exponentials appearing in the RHS of (15), since this is the opposite of that of the corresponding functions contained in the RHS of (12).
The derivation of the CSDE follows the same line of reasoning as the CSFE. Moreover, the latter algorithm has the same structure as the former one, because of the similarity between the involved signal models. The few differences, originating from the above mentioned sign reversal, can be summarised as follows: 1) The vectors X 2) The estimateˆ (i) of the delay residual is evaluated on the basis of (121) and (119), where, however, andÂ is the estimate of the complex amplitude A computed on the basis of (125) (the coefficient c(ρ), instead, is still expressed by (127)). Moreover, in the evaluation of the coefficients {b(ρ), c(ρ)}, the quantitȳ X k,ρ (114) is replaced bȳ whereρ is still defined by (122). Given the CSDE, the counterpart of the CSFEC algorithm, called complex single delay estimation and cancellation (CSDEC), can be easily developed.

D. TARGET CANCELLATION PROCEDURES EMPLOYED IN FMCW RADAR SYSTEMS
In T2 of the RASCA-FR2 and RASCA-FR3 (and, in particular, in STDREC-S2; see (53)), a target cancellation procedure is used in combination with the SFE. This procedure requires the evaluation of the triad (C (v) , that represents the contribution given by the i-th (i.e., by the last) point target detected on the v-th VA. IfF denote the estimates of the normalised frequency and the complex amplitude, respectively, characterizing this target, the expressions and with k = 0, 1 and 2,W is the normalised frequency associated with the frequencyf It is important to point out that an efficient method can be used for the computation of the vectorsW respectively, where Therefore, the identities listed in [10, eqs. (84)- (85) and (145)] can be exploited for an efficient computation of the RHSs of (136) and (137). A target cancellation procedure is also employed in T2 of the RASCA-FC2 and RASCA-FC3; however, in this case, the CSFE is adopted in place of the SFE, and the vectors C and with k = 0, 1 and 2, andw (v) i is still expressed by (134). The vectorW (v) k [i] appearing in (130)-(132) (with k = 0, 1 and 2) can be efficiently computed following the same approach illustrated above for the SFE.
The CSFE is also employed in T3-S1 and, in particular, in STDAEC-S2 of the RASCA-FR2, RASCA-FR3, RASCA-FC2 and RASCA-FC3. In this case, the cancellation procedure requires the evaluation of the contribution given by the i-th (i.e., by the last) target detected in the l-th frequency bin to the whole array (see (62)). Here, we focus on the target cancellation procedure employed in the above mentioned RASCAs. In this case, denote the estimates of the complex amplitude, the normalised vertical spatial frequency and the normalised horizontal spatial frequency, respectively, characterizing the i-th target, the expression is employed for any VA (i.e., for any p and q). Finally, it is important to mention that the cancellation procedure adopted in STDREC algorithm aims at removing the contribution of a single target in each of its iterations. If a cluster of m (v) i distinct frequencies is estimated by the SFE (CSFE) in the i-th iteration of the above mentioned algorithm, each of the components of the triad (C (v) i terms and each term is evaluated on the basis of (130)-(132) ((140)-(142)).

E. TARGET CANCELLATION PROCEDURES EMPLOYED IN SFCW RADAR SYSTEMS
In this paragraph we take into consideration a MIMO SFCW radar system and focus on the two cancellation procedures employed in the RASCA-S2 and RASCA-S3 (see Paragraph V); the first procedure is used in combination with the CSDE in T2 (and, in particular, in STDREC-S2), whereas the second one with the CSDE in T3-S1 (and, in particular, in STDAEC-S2). The first procedure is based on (53); in the considered radar system, the vectors (C (v) appearing in the RHS of that equation and expressing the contribution given by the i-th (i.e., by the last) target detected on the v-th VA are evaluated as (147) VOLUME 10, 2022 and C (v) with k = 0, 1 and 2, represent the estimates of the complex amplitude A whereq Even in this case, the identities listed in [10, eqs. (84), (85) and (145)] can be exploited for an efficient computation of the RHS of (152). The second cancellation procedure is expressed by (62) and requires the evaluation of the matrix C for any VA (i.e., for any p and q).

VII. LIMITATIONS
In this section, some technical limitations that have emerged in the implementation of our algorithms on commercial radar devices are illustrated and the solutions we have devised to mitigate their impact are described.

A. UNEQUAL RESPONSE OF VIRTUAL ANTENNAS
The derivation of the RASCAs for FMCW radar systems relies on the assumption that the real (complex) sample sequence made available by the v-th VA is expressed by (6) ((12)); the validity of a similar model has been assumed for SFCW radar systems (see (15)). The adopted signal models hold if the amplitudes of the L overlapped oscillations contributing to the useful component of the received signal do not change from antenna to antenna. However, our experiments accomplished on commercial colocated radar devices have evidenced that: a) these amplitudes are not constant across the whole virtual array; b) their differences are influenced by the azimuth and the elevation of each target. We believe that all this is due to the different behavior of the multiple receive chains employed in each MIMO device and to the mismatches in the receive antenna patterns. This problem affects the quality of the data acquired through both FMCW and SFCW radar systems and, consequently, the accuracy of our detection and estimation algorithms. It can be mitigated by enriching the physical array with a set of surrounding passive antennas; in this case, the array is artificially extended with new antennas along all its sides, so that the behavior of all its active antennas becomes more uniform.
It important to point out that, in principle, the presence of this phenomenon can be accounted for in the development of target detection and estimation algorithms by including its effects in the received signal model. For instance, (6) can be generalised as where α v (θ l , φ l ) represents an attenuation factor depending on the angular coordinates of the l-th target and v is the VA index. Consequently, the complex amplitude associated with the l-th target detectable on the considered VA becomes (see (7)) Neglecting the presence of the factor α v θ l , φ l in the development of our algorithms has the following implication: an error is introduced by the STDAEC algorithm in its cancellation procedure (see STDAEC-S2 in Paragraph VI). Note, in particular, that the estimateĈ i [l] of the complex amplitude characterizing to the i-th target detected in theα l -th bin is computed after the overall spatial folding (i.e., after STDAE-S5); consequently, its absolute value represents a sort of spatial average computed over all the involved VAs. Moreover, only the phase variations of this complex gain are accounted for in the computation of the contribution C in the evaluation of the term C (i) X 0 [l] appearing in (144)-(145). Estimating the function α v (θ, φ), however, is a time consuming task, since it requires a proper measurement setup and an anechoic chamber. We believe that this problem can be circumvented by: a) exploiting deep learning techniques [27] in the SPE; b) adopting a data-driven approach [28], [29]. This solution is motivated by the fact that: a) Deep learning techniques can be employed to approximate complicated functions, that do not lend themselves to a simple parametric representation and without requiring particular expertise in data pre-processing. b) A data-driven approach allows to train different models on the basis of data collected in a real scenario or synthetically generated data, without prior knowledge about the parametric representation of the considered problem. Note that a fundamental role is played by the adopted training procedure since it makes the involved network able to generate correct predictions on the basis of never seen data available at its input.
In practice, the adoption of the proposed approach requires modifying the STDAEC technique employed in the RASCAs (see Fig. 5) and, in particular, embedding a deep neural network in it. This network is employed to estimate the distorted amplitudes of all the targets detected in the l -th frequency bin (with l = 0, 1, . . . , L b − 1), so that accurate cancellation becomes possible.
The use of this solution in our radar systems is not investigated in the following, since it is out of the scope of this manuscript.

B. ANTENNA COUPLING
In our description of the SFE, the CSFE and the CSDE (see Section VI), it has been implicitly assumed that the minimum frequency of the useful component contained in the observed data sequence can be arbitrarily small. Unluckily, this is not always true. For instance, in commercial colocated FMCW MIMO radar systems, a strong interference is observed in the lower portion of the spectrum evaluated on all the receive antennas. This phenomenon, known as mutual coupling [30], is due to the electromagnetic coupling that originates from the small distance between adjacent transmit and receive antennas [20]. Its impact can be mitigated resorting to various methods based on calibration measurements [31]. Because of mutual coupling, any target whose range is below a certain threshold cannot be detected by our algorithms in a reliable fashion.

VIII. OTHER TARGET DETECTION AND ESTIMATION TECHNIQUES
The detection and estimation algorithms described above have been compared, in terms of accuracy and complexity, with two different types of algorithms that, similarly as the RASCAs, are able to generate radar images in the form of point clouds. The algorithms of the first type are called FFT-based algorithms (FFT-BAs), since they rely on multidimensional FFT processing for the evaluation of all the spatial coordinates of targets (i.e., their range and DOA); such algorithms have been inspired by the FFT-based algorithm proposed by Texas Instrument in [5]. The algorithms of the second type, instead, are called MUSIC-based algorithms (MUSIC-BAs); these make use of the same method as the first type for range estimation, but the MUSIC algorithm for DOA estimation [32]. In the remaining part of this section, a brief description is provided for both types. In our analysis, we always refer to a FMCW radar system, since, if a SFCW radar system is considered, the only change to be made in the described algorithms consists in replacing each FFT with an IFFT of the same order.
The inner structure of both types of algorithms is described by the block diagram shown in Fig. 9. The processing accomplished by the blocks this diagram consists of, can be summarized as follows. Each vector of the set {x  (48)). Based on this set of vectors, the N 0 -dimensional power spectrum is computed; here, with i = 0, 1, . . ., N 0 − 1. The vector P 0 (158) feeds the cell-averaging smallest of -constant false alarm rate (CFAR-CASO) algorithm developed in [33]. Based on this algorithm, a target is detected in the i-th frequency bin if where i ∈ {i m , i m + 1, . . ., i M }. Here, represents a decision threshold, K 0 is a real parameter whose value is selected on the basis of the required false alarm rate, VOLUME 10, 2022 represent the average of the power spectrum computed over C s adjacent bins positioned on the left and on the right, respectively, with respect to the i-th frequency bin. Moreover, G s and C s are two integer parameters defining the size and the position (with respect to the i-th bin), respectively, of the set of frequency bins involved in the computation ofP l (162) and P u (163 ), whereas i m and i M are two non negative integers such that i m ≥ (G s + C s ) and i M ≤ N 0 − 1 − (G s + C s ). In our work, the inequality is also required to be satisfied together with the condition (160), where P l,u represents the largest element of the This allows us to reduce the overall number of detected targets, so reducing the density of the generated point cloud.
The CFAR-CASO algorithm generates the vector whereα l represents the index of the frequency bin in which the l-th target has been detected (with l = 0, 1, . . . , L b − 1) and L b is the overall number of detected targets. This vector is processed for DOA estimation. The two options (associated with the above mentioned types of algorithms) are considered for this task and are described in the remaining part of this paragraph.
FFT-based DOA estimation -Let us focus first on the case in which a virtual HULA, consisting of N VH virtual elements, is employed for resolving the targets associated with a given frequency bin and estimating their azimuth. In this case, azimuth estimation consists of the following two steps: 1) The N VH -dimensional column vector (see (88)) collecting the spectral information available on the whole array and referring to theα l -th frequency bin (with l = 0, 1, . . ., L b − 1) is applied to anN 0 order FFT algorithm; let 2) The dominant peaks 11 in the sequence {|s k [l]|; k = 0, 1, . . . ,N 0 − 1} are identified; each peak corresponds to a 11 It is important to distinguish peaks associated with different targets from side-lobes; in our simulations, a candidate peak is classified as a side-lobe (and, consequently, ignored) if its amplitude differs by more than 1 dB from that of a close dominant peak, as suggested in [5]. distinct target. If k i [l] denotes the index of i-th peak (with i = 0, 1, . . ., L h [l] − 1, where L h [l] is the overall number of targets detected in the considered frequency bin), the estimate of the azimuth of the i-th target is evaluated aŝ where Let us now consider the case in which the URA represented in Fig. 4 is employed for resolving the targets associated with each frequency bin, and estimating their azimuth and elevation. The algorithm employed in this case involves the [X 0,α l [p, q]] (60), collecting the spectral information available on the whole array for theα l -th frequency bin. This algorithm consists of the following four steps: 1) The p R -th row of the matrix X[l] is processed to generate theN 0 -dimensional column vector s VULA,0 [l] = [s 0,0 [l] , s 0,1 [l] , . . . , s 0,N 0 −1 [l]] T on the basis of (68); here, p R represents the column index of the reference antenna in the considered URA (see Fig. 4).
2) The dominant peaks of the sequence is the overall number of targets detected in the considered frequency bin), the estimate of the elevationφ i [l] of the associated target is evaluated asφ 3) The 2D FFT of the matrix X[l] is computed; this produces theN 0 ×N 0 matrixS[l] = [S k,r [l]], such that where  The last step concludes our description of the FFT-BAs. Note that the overall number of detected targets is given bŷ MUSIC-based DOA estimation -Similarly as our description of the FFT-BAs, we first focus on the case in which a virtual HULA, consisting of N VH virtual elements, is employed for resolving the targets associated with a given frequency bin and estimating their azimuth. In this case, the algorithm considered for DOA estimation consists of the following three steps: 1) The N VH × N VH autocorrelation matrix is computed; here, X [l] is defined by (166).
2) TheN 0 -dimensional pseudo-spectrum P (l) MU is evaluated; its k-th element is given by with k = 0, 1, . . . ,N 0 − 1; here, (·) H denotes the conjugate and transpose operator, Q N VH is a matrix having size N VH × (N VH −1) and whose columns are the (N VH −1) noise eigenvectors (associated with the (N VH −1) smallest eigenvalues) of R X [l] (174) and a[k] is a N VH -dimensional steering vector, whose n-th element a n [k] is given by 12 a n [k] = exp jπ n hN with n = 0, 1, . . . , N VH − 1.
3) The dominant peaks appearing in the sequence {P Let us consider now the case in which the uniform rectangular array shown in Fig. 4 is employed for resolving the targets associated with each frequency bin, and estimating their azimuth and elevation. In this case, the adopted procedure involves the N VH × N VV matrix X[l] [X 0,α l [p, q]] (60) for anyα l and consists of the following four steps: 1) The pseudo-spectrum referring to the reference VULA (that consists of N VULA virtual elements) is evaluated. In this step, we assume that the p R -th row of X[l] is employed for the evaluation of the autocorrelation matrix R X [l] (174) and that theN 0 -dimensional vector P (VULA) MU [l] is computed on the basis of (175)-(176) (note that N VR and δ[k] are replaced by N VULA and δ[r], respectively).
2) The dominant peaks appearing in the sequence of the ordered elements of P (VULA) MU [l] are identified; let L v [l] denote 12 Note that the following equation applies to an FMCW radar; if an SFCW radar is considered, the sign of the argument of the complex exponential appearing in its RHS must be reversed. their overall number. If the i-th peak is found for r = r i [l] (with i = 0, 1, . . ., L v [l] − 1), the elevationφ i [l] of the associated target is evaluated on the basis of (169).
3) The pseudo-spectrum P (HULA) MU [l, i] associated with the ith estimated elevation is evaluated for the whole virtual array. In this step, if we assume that the autocorrelation matrix R X is computed according to (174) (where, however, X [l] is the N VH ×N VV matrix defined above), theN 0 -dimensional vector P (HULA) MU [l, i] is generated on the basis of (175). Note that, in this case, N VR is replaced by N HULA and that the n-th element a n [k] of the N HULA -dimensional steering vector a[k] is The order adopted in the computation of the pseudo-spectra (first the vertical pseudo spectrum P The performance of the FFT-BAs and the MUSIC-BAs has been assessed for both FMCW and SFCW radar systems working in 2D and 3D propagation scenarios. The acronyms adopted in the following for these types of algorithms are summarized in Table 2.

IX. COMPUTATIONAL COMPLEXITY
The computational cost of the algorithms described in Sections IV, V and VIII has been carefully assessed in terms of floating point operations (flops) to be executed in the detection of L targets. 13 Various details about the method we adopted for the evaluation of the computational cost of each algorithm are provided in Appendix for the RASCA-FC3 only. Our analysis leads to the conclusion that the overall cost of this algorithm and RASCA-FC2 is approximately of order O(M R−FC3 ) and O(M R−FC2 ), respectively, where (see (195)) and and M M−FC2 = 8N VH N 0 log 2 (N 0 )+L bN0 (N 3 VH +16N 2 VH ). (184) It is important to keep in mind that a comparison among the computational costs listed above does not fully account for the gap that can be observed in the execution speed of the corresponding algorithms. In fact, in practice, a portion of the computation time is absorbed by the procedure employed to find the dominant peaks of real sequences in both the FFT-BAs and the MUSIC-BAs. Moreover, the vector A CF (165), collecting the indices of the frequency bins in which at least one target has been detected, may include ghost targets; as evidenced by our computer simulations, the impact of this phenomenon on the overall computation time may not be negligible. Despite this, some interesting insights on how the complexity is influenced by the overall number of targets can be obtained by comparing the computational costs (179), (181) and (183) ((180), (182) and (184 )) in two specific scenarios. The first scenario we take into consideration refers to the case in which the mutual distance between the targets is above the range resolution of the employed radar system, so that K T 2 = L, K T 3 = 1 and L b = L can be assumed in the RHS of (179)-(184). In our second scenario, instead, the targets form clusters, each of which consists of four targets having the same range, but different angular coordinates; for this reason, K T 2 = L/4, K T 3 = 4 and L b = L/4 can be assumed in the RHS of (179)-(184). Moreover, the following parameters have been chosen for both scenarios: a) N VR = 256; b) N 0 = 1024; c) N A = 10 d) N VV = 16; e) N VH = 16; f)N 0 = 32. The dependence of the complexity M alg on L is represented in Fig. 10a (Fig. 10b) for the first (second) scenario; here, alg denotes the algorithm which this complexity refers to. From these figures it is easily inferred that: a

X. NUMERICAL RESULTS
In this section, the accuracy of the RASCAs is assessed on the basis of both synthetically generated and experimental data, and is compared with that provided by various FFT-BAs and MUSIC-BAs.

A. NUMERICAL RESULTS BASED ON SYNTHETICALLY GENERATED MEASUREMENTS
In this paragraph, the accuracy achieved by the RASCA-FC3, the FFT-FC3 and the MUSIC-FC3 in the generation TABLE 3. Root mean square errorε X , peak errorε X and detection rate R D evaluated in the two simulation scenarios defined in Paragraph X-A. Target range, azimuth and elevation are taken into consideration. Note that, in principle, the available antenna array allows us to achieve the range resolution the azimuthal resolution and the elevation resolution The considered radar system is assumed to operate in the presence of L = 10 targets, whose echoes have unit amplitude. The range, the azimuth and the elevation of each target are sequentially generated at the beginning of each run. Moreover, the range R k , the azimuth θ k and the elevation φ k of the k-th target (with k = 1, 2, . . . , 10) have been randomly evaluated in a way that: a) they belong to the intervals [1,10] m, [−π/3, π/3] rad and [−π/3, π/3] rad, respectively; b) the minimum spacing between the k-th target and the previously generated (k − 1) targets is not smaller than R (185), θ (186) and φ (187) in the range, azimuth and elevation dimensions, respectively (scenario S1) or is not smaller than R (185) in the range domain, but can be arbitrarily small in the azimuth and elevation dimensions (this scenario is denoted S2). In our computer simulations, here, X i andX i represent the exact value of a parameter X and its corresponding estimate, whereas N m represents the overall number of synthetically generated values of X ; note that, if all the targets are detected by the considered algorithm in each run, 190) where N r is the overall number of simulation runs. In our work, the performance of the above mentioned algorithms has been assessed by: a) evaluating the detection rate for both the considered scenarios; b) ignoring the failure events (i.e., the events in which not all the targets have been detected) in the evaluation of all the RMSEs. The three performance indices defined above have been assessed on the basis of the estimates generated by executing N r = 500 runs; the resulting values are summarised in Table 3. From these results it is easily inferred that: a) The RASCA-FC3 achieves the lowest RMSEs in range, azimuth and elevation (range, elevation) estimation in the first (second) scenario; for instance, the RMSEε θ characterising the RASCA-FC3 is about 1.3 (1.3) times smaller than the corresponding RMSE obtained for the FFT-FC3 in the first (second) scenario. b) The RASCA-FC3 exhibits the lowest peak errors in range, azimuth and elevation (range) in the first (second) scenario; for instance, its peak error ε R is 2 times smaller than the corresponding RMSE obtained for the FFT-FC3 and the MUSIC-FC3 in both scenarios. c) All the considered algorithms achieve an excellent accuracy in both scenarios, since the RMSEs evaluated for range, azimuth and elevation are smaller than the corresponding resolutions given above. d) The FFT-FC3 and the MUSIC-FC3 are outperformed by the RASCA-FC3 in terms of detection rate; in fact, the value of this parameter is about 70 % for the first two algorithms, but is equal to 100 % for the RASCA-FC3, since the last algorithm has been able to detect all the targets in every simulation run in both scenarios.

B. NUMERICAL RESULTS BASED ON EXPERIMENTAL MEASUREMENTS
In this paragraph, we first describe the radar devices employed in our measurement campaigns and the adopted experimental setup. Then, we analyse: 1) the accuracy achieved by our RPE (and, in particular, by the STDREC algorithm) in range and phase estimation on multiple antennas of the same array in the presence of a single target and of multiple targets; 2) the accuracy of the 2D (3D) images generated by RASCA-FR2 (RASCA-FR3), RASCA-FC2 (RASCA-FC3) and RASCA-S2 (RASCA-S3) in the presence of multiple targets.

1) EMPLOYED RADAR DEVICES AND ADOPTED EXPERIMENTAL SETUP
A measurement campaign has been accomplished to acquire a data set through two FMCW MIMO radars and one SFCW MIMO radar, all operating in the E-band. The first FMCW device, dubbed TI FMCW radar in the following, is the TIDEP-01012 Cascade mmWave radar (see Fig. 11 -a)). It is manufactured by Texas Instrument Inc. [34], classified as a long range radar (LLR) and provides both the in-phase and quadrature components of received signals (i.e., complex measurements). Its main parameters are: a) chirp slope µ = 4 · 10 13 Hz/s; b) bandwidth B 1 = 2.5 GHz; c) central frequency f c = 77 GHz; d) sampling frequency f s = 8 MHz; e) number of samples per chirp N = 512. Moreover, it is endowed with a planar array made of N T = 12 TX and N R = 16 RX antennas (each consisting of an array of four patch elements), as shown in Fig. 11-a). The corresponding virtual array consists of 12 · 16 = 192 VAs; however, only 134 of them are available, since the remaining 58 VAs overlap with the other elements of the virtual array. As shown in Fig. 11-b) (where each VA is represented by a small blue circle), the virtual array has the following characteristics: 1) the non-overlapped VAs form an horizontal ULA (HULA 1 ), consisting of N HULA 1 = 86 VAs and three smaller HULAs, each made of 16 equally-spaced VAs; 2) the inter-antenna spacing of all the HULAs is d VH = λ/4; 3) the vertical spacing of the three smaller HULAs is not uniform, since d VV 1 = λ/4, d VV 2 = λ and d VV 3 = 3λ/2 (see Fig. 11-b)). This virtual antenna array allows us to achieve a range, azimuth and elevation resolution equal to R 1 = 5.8 cm (see (185)), θ 1 = 1.35 • (see (186)), and φ 1 = 16.4 • (see (187)), respectively; note that the elevation resolution is coarser than the azimuth one since N VV = 7 equally aligned antennas (d VV = d VV 1 = λ/4) are assumed along the vertical direction (this is equivalent to considering an elevation aperture D y = 3λ along the vertical direction; see [20]). In our work, on the one hand, a central portion of the first HULA (contained inside the red rectangle appearing in Fig. 11-b)), consisting of N VV = 16 antennas, has been exploited for 2D imaging, in order to guarantee a fair comparison with the other two radar devices. On the other hand, the whole array have been employed for 3D imaging.
The second FMCW device, dubbed Inras FMCW radar in the following (see Fig. 12-a)), is a modular system manufactured by Inras GmbH [35] and consisting of: a) the so called Radar Log board; b) an RF front-end including multiple TX/RX antennas and monolithic microwave integrated circuits (MMIC) operating at 77 GHz. This system is classified as a LLR and its main parameters are: a) chirp slope µ = 9.7656 · 10 12 Hz/s; b) bandwidth B 2 = 2.5 GHz; c) central frequency f c = 77 GHz; d) sampling frequency f s = 8 MHz; e) number of samples per chirp N = 2048. Unlike the TI FMCW radar, this device provides only the in-phase component of the RF received signals and, consequently, real measurements. Moreover, it is endowed with a custom designed planar array made of N T = 16 TX antennas and N R = 16 RX antennas, each consisting of an array of six patch elements, as shown in Fig. 12-a). The resulting virtual array, consisting of N VR = 16 · 16 = 256 VAs is shown in Fig. 12-b). As it can be inferred from the last figure, the virtual array has the following characteristics: 1) It consists of 16 HULAs, each of which is made of 16 antennas with inter-antenna spacing d VH = λ/4. 2) The vertical distance between each couple of its adjacent HULAs is d VV = λ/2; this entailes the unambiguous elevation range [−45 • , 45 • ]. 3) Its shape is not rectangular (the horizontal shift of adjacent HULAs is equal to λ/4). This virtual array allows us to achieve the same range resolution as the TI FMCW radar, and azimuth and elevation resolutions equal to θ 2 = 7.6 • and φ 2 = 3.8 • , respectively (see ( 186)-(187)). In our work, the HULA contained inside the red rectangle appearing in Fig. 12-b) (the whole array) has been exploited for 2D (3D) imaging.
The third device, dubbed VIC SFCW radar in the following, is the Vito-In-Car radar manufactured by Vayyar Imaging Ltd Company [36] (see Fig. 13-a)); it is classified as a short range radar (SRR). Its main parameters are: a) initial transmit frequency f 0 = 78 GHz; b) frequency spacing f = 16.67 MHz; c) overall number of frequencies N f = 121. Therefore, the bandwidth and the central frequency of the radiated signal are B 3 = 2.0 GHz and f c = 79 GHz, respectively. This device is equipped with the URA shown in Fig. 13-a); this is composed by N T = 16 TX antennas and N R = 21 RX antennas, so that an URA of 16 · 21 = 336 VAs is available, as shown in Fig. 13-b). However, in our work, only N VV = 20 HULAs, each consisting of N HULA = 16 virtual channels characterized by the inter-antenna spacing d VH = λ/4, have been exploited; note that the vertical distance between two adjacent HULAs is d VV = λ/4. The available array, made of N VV = 320 virtual elements, allows us to achieve a range resolution slightly better that that provided by the TI FMCW radar ( R 3 = 7.5 cm; see (185)), and azimuth and elevation resolutions equal to θ 3 = 7.6 • and φ 3 = 6.0 • , respectively (see (186)-(187)). In our work, the HULA contained inside the red rectangle appearing in Fig. 13-b) (the whole array) has been exploited for 2D (3D) imaging.
Our measurement campaigns have been conducted in a large empty room (whose width, depth and height are 10 m, 8 m and 2.5 m, respectively). Each of the employed radar devices has been mounted on an horizontal wooden bar together with a pico-flexx camera manufactured by PMD Technologies Inc. [37] and has been lifted by a tripod at an height of roughly 1.60 m from ground, as shown in Fig. 14. The employed camera is based on a near-infrared vertical cavity surface emitting laser, and is able to provide a depth map or, equivalently, a 3D point-cloud of a small region of the observed environment (its maximum depth is equal to 4 m, whereas its field of view is 62 • × 45 • ). In each measurement campaign, the experiments have been repeated for all the radar devices exactly in the same conditions. Note that the measurements acquired through the SFCW MIMO radar have been pre-processed by an on-board cancellation algorithm that exploits the measurements acquired from the first transmitted frame to cancel out unwanted received echoes.
It is important to point out that: a) in all the radar systems, the target ranges have been estimated with respect to the central virtual channel of the employed ULA; b) the exact target positions have been acquired with respect to the centre of the pico-flexx camera; c) the data processing has been accomplished in the MATLAB environment; d) all our detection and estimation algorithms have been run on a desktop computer equipped with a single i7 processor.

2) RANGE AND AMPLITUDE ESTIMATION
In this paragraph, the accuracy of the STDREC algorithm employed by the RPE is analysed for two specific static scenarios. The first scenario is characterized by a single detectable target (a small metal disk 14 having a diameter equal to 5.5 cm) placed in ten different positions. The target range R and azimuth θ have been selected in the interval  Table 4 for all the employed radar devices (the data referring to the i-th position are collected in the column identified by T i , with i = 1, 2, . . ., 10). The second scenario, instead, is characterized by the presence of an overall number of targets ranging from 1 to 9 (so that 1 ≤ L ≤ 9). The targets are represented by 14 Each target is hung from the ceiling: a nylon thread has been used for suspending it. small coins with a diameter of 2 cm; the range and azimuth characterizing their exact positions are listed in Table 5 (the data referring to the i-th target are collected in the column identified by T i , with i = 1, 2, . . ., 9). Each target has been sequentially added in our scenario; this has allowed us to assess how the performance of the STDREC algorithm is influenced by the value of the parameter L in the presence of closely spaced targets.
Prior knowledge of L has been assumed during the processing; moreover, the following values have been selected for the parameters of the STDREC algorithm 15  Note that: a) the value of the oversampling factor (M ) has been selected in way to guarantee approximately the same value of N 0 in all cases, i.e. roughly the same resolution in the spectral analysis of radar signals; b) the values of the parameters N SFE , N CSFE and N CSDE are all equal and large enough so that accurate range estimation is achieved by the STDREC algorithm.
The accuracy of range estimates has been assessed by evaluating the RMSEε R and the peak errorε R , expressed by (188)-(189) with X = R, X i = R i andX i . Since the RCS of the considered targets was unknown, our analysis of the complex gains available over the 16 channels of the considered virtual ULA and associated with the same target has concerned only their (unwrapped) phase. The phases {ψ (v) ; v = 1, 2, . . . , 16} estimated by the STDREC algorithm over the considered reference HULA (consisting of 16 VAs; see the red rounded rectangles appearing in Figs. 11-b), 12-b), 13-b)) and associated with a target placed at approximately 16 the 15 Note that in this case the stopping criterion based on eq. (54) has not be employed, since the overall number of targets is known. 16 The exact range of this target can be found in the T 7 (T 3 ) column for the TI FMCW and the VIC SFCW radars (Inras FMCW radar) in Table 4.  same azimuth angle with respect to the centre of the radars is shown in Fig. 15. Since the distance d VH between adjacent virtual channels is constant, the (unwrapped) estimated phases exhibit a linear dependence on the index of the virtual channel, as illustrated in Section II (see, in particular, (10) and (11)). Moreover, if a linear fitting is drawn for these data, it should be expected that the slope of the resulting straight line is proportional to sin(θ) (see (23) with φ = 0); this is confirmed by the results shown in Fig. 15 for both the FMCW and the SFCW radar systems. To assess the quality of the estimated phases, their RMSEε ψ has been evaluated in all the scenarios; in doing so, the linear fitting of the 16 phases } has been taken as a reference with respect to which the error of each of them has been computed.
The estimate of the target range generated by the STDREC algorithm for each of the N m = 10 distinct positions considered in the first scenario are listed in Table 4; in the same table, the value ofε ψ computed for each position is also given. The target ranges and their estimates listed in Table 4 are also represented in Fig. 16. The errorsε R andε R , the mean ofε ψ (denotedε m,ψ and generated by taking the average of the N m values available forε ψ ) and the average computation time (CT) evaluated on the basis of these results are listed in Table 6.
The results referring to the first scenario lead us to the following conclusions: 1) In all the considered cases, the STDREC is able to accurately estimate the range and the phase characterizing each target. 2) All the values ofε R andε R are comparable, reasonably low and in the order of the resolution of our devices.
3) The range estimates evaluated for the VIC radar are closer to the corresponding exact values when the target range is below 3 m, while the opposite is true for the other two devices; this can be related to the fact that the first device is a SRR, whereas the other two are LRRs 4) The Inras FMCW radar achieves the lowestε m,ψ . 5) The value ofε ψ characterizing the VIC SFWC radar is lower when the target is closer to this device, whereas the opposite occurs in the case of the TI FMCW radar. This result is mainly due to: a) the far-field condition is satisfied for the VIC SFWC device at shorter ranges than for the TI FMCW device, since the antenna array of the former is more compact than that of the latter; b) the VIC SFCW (TI FMCW) radar is a SRR (LRR) system. For these reasons, the measurements provided by the VIC SFCW radar allow us to achieve more accurate results at shorter distances. On the contrary, the higher transmit power radiated by the TI FMCW radar explains why a smallerε ψ is achieved at larger ranges. 6) The CTs are always in the order of few milliseconds. Let us focus on the second scenario. In this case, our range estimates have been generated by: a) the STDREC algorithm for all the radar devices; b) the HSTDREC algorithm for the TI FMCW and the VIC SFWC radars; c) the Alg-YA, the Alg-S and the CFH algorithm for the TI FMCW. The obtained results are listed in Table 5. The errorsε R andε R , and the CT obtained in this case are listed in Table 7. From these results it can be inferred that: 1) In the case of the TI FMCW radar, all the considered algorithms achieve comparable accuracy. However, the STDREC and the HSTDREC algorithms, unlike all the other algorithms, achieve the lowest RMSE and peak error.
2) The HSTDREC algorithm is not more accurate than the STDREC algorithm; moreover, these algorithms are characterized by similar CTs.
3) The estimated RMSEs and peak errors are in the order of the resolution of our radar devices, but a little bit higher in the Inras FMCW radar systems. This is mainly due to the poorer estimates evaluated for the targets T 8 and T 9 , since, in our specific experiment, the energy received from these targets has been found to be lower than that coming from the others. This problem is  not so evident in the case of the TI FMCW and SFCW radars, whose RMSEs and peak errors are very low. Finally, we would like to stress that the accuracy of STDREC and HSTDREC algorithms can be related to the accuracy of the estimation and cancellation procedure they accomplish. This is exemplified by Fig. 17, where the initial amplitude spectrum of the signal received on the central virtual channel of the TI FMCW radar in the second scenario and its (weak) residual, resulting from the cancellation of the spectral contributions due to the detected targets, are shown. Here, the range and amplitude of the targets estimated by the STDREC (HSTDREC) are also represented by red circles (green crosses).

3) TWO-DIMENSIONAL AND THREE-DIMENSIONAL IMAGING
In this paragraph, the accuracy of the 2D and 3D images generated by the RASCAs is assessed. Two different groups of experiments have been carried out. The first (second) group of experiments has allowed us to assess the performance achieved by the above mentioned algorithm in 2D (3D) imaging. In both cases, the measurements have been acquired in the presence of an increasing number of targets for all our radar devices. In the first group of experiments, the following choices have been made:    Table 8).
2) The measurements have been acquired through a virtual ULA, consisting of 16 VAs, in all the considered radar systems. As far as the second group of experiments is concerned, the following choices have been made: 1) The range, azimuth and elevation of the targets have been selected in the intervals [1.9, 2.8] m, [−30 • , 35 • ] and [−10 • , 10 • ], respectively (see Table 11).
2) The measurements have been acquired through the whole virtual array of each of our radar devices. The following values have been selected for the parameters of the RASCAs: a) N A = 16 (N A = 10) in the RPE employed in 2D (3D) imaging; b) N CSFE = N SFE = N CSDE = 5 in both the STDREC and the STDAEC algorithms; c) the threshold T OF = 0 has been selected in the STDAE-S4 algorithm; d) the values of the parameters N 0 and M are equal to those employed for the STDREC in the previous paragraph; e) the oversampling factor isM = 16 (M = 7) for Inras FMCW and VIC SFCW radar (for the TI FMCW), so that the FFT orders areN 0 = 256 andN 0 = 320 (N 0 = 602), respectively. Moreover, the following values have been selected for the parameters 17 of the FFT-BAs and the MUSIC-BAs: C s = 6, G s = 6 and K 0 = 2. Prior knowledge of L has been assumed and the threshold T STDAEC has been selected in the range [0.01, 0.9] · E (0) [l] (63) (the value of this threshold has been adjusted on the basis of the SNR characterizing the received signal and the overall number of detectable targets).
The estimates of range and azimuth generated by the RASCAs on the basis of the measurements acquired in our first group of experiments are listed in Table 8, whereas the values of RMSE, peak error and CT computed by averaging the RMSEs, peak errors and CTs evaluated in each single experiment are listed in Table 9. In the last table, the values of RMSE, peak error and CT for the employed FFT-BAs and MUSIC-BAs are also provided. These results lead to the following conclusions: 1) All the range and azimuth errors are comparable with the resolution of our devices. 2) The RASCAs always outperform the other algorithms and require a lower CT.
3) The highest range (azimuth) peak errors and RMSEs are found in the case of the Inras FMCW (VIC SFWC) radar; the TI FMCW radar, instead, achieves the lowest range and azimuth errors. This is mainly due to the differences in the SNR available at the receive side of distinct radar devices is different. The good accuracy achieved by the RASCAs is also evidenced by Fig. 18, where a range-azimuth map [19], generated though standard 2D FFT processing of the measurements acquired through the Inras FMWC radar, is represented as a contour plot 18 ; in the same figure, the exact position of the five targets employed in our first group of experiments and their estimates evaluated by all the considered algorithms are shown.
Let us consider now on the results obtained for our second group of experiments. The estimates of range, azimuth and elevation generated by the RASCAs are listed in Table 11, whereas the values of RMSE, peak error and CT evaluated on the basis of this table are listed in Table 10. In the last table, the errors characterizing the FFT-based and MUSIC-based algorithm for 3D imaging are also provided. From these results it can be inferred that: 1) The RMSEs and the peak errors evaluated for target range, azimuth and elevation are reasonably low and comparable with those obtained in the case of 2D imaging. Moreover, these errors are smaller than the ones characterizing the FFT-based and MUSIC-based algorithms. 17 Our simulations have evidenced that small changes in the value of these parameters do not significantly influence the detection probability and the estimation accuracy of the considered algorithms. 18 Note that x − y coordinates are employed in this case, in place of range and azimuth; the position of the radar system corresponds to the origin of our reference system. VOLUME 10, 2022 TABLE 8. Exact range and azimuth of the five targets considered in our first group of experiments and corresponding estimates generated by the RASCAs.

TABLE 9.
Root mean square errorε X , peak errorε X , and computation time (CT) evaluated on the basis of our first group of measurements. Target range and azimuth are taken into consideration.
2) The range errors evaluated for all the radar devices are comparable; however, the highest (lowest) peak error  Table 8) and the other algorithms are also shown. The rectangles allow to delimit the region in which the position of each target and its estimates are located. and RMSE are achieved in the case of the VIC SFCW (TI FMCW) radar.
3) The azimuth and elevation estimates computed on the basis of the measurements acquired through the TI and Inras FMCW radars are reasonably good. 4) The accuracy achieved by the VIC SFCW radar in the estimation of target elevation is better than that in azimuth estimation. This is due to the fact that the URA of this device provides a larger number of channels (i.e., a finer angular resolution) along the vertical direction with respect to the horizontal one. 5) The average CT is in the order of few seconds for all the proposed algorithms; the lowest (highest) average CT is found in the case of the VIC SFCW and the TI FMCW radars (Inras FMCW radar). This is mainly due to the fact that the STDREC algorithm employed in the case of the Inras FMCW requires an higher computational effort with respect to its counterpart employed with the TI FMCW and the VIC SFCW radars. 6) In general, the CT of RASCAs is higher than that required by the proposed FFT-based or MUSIC-based methods. This is mainly due to the recursive cancellation procedure that is not employed by the other two methods. However, we believe that this cancellation procedure plays a fundamental role in the detection of weak targets and allows to achieve a good estimation accuracy. The good accuracy and resolution provided by the RAS-CAs are highlighted by Fig. 19, where the exact positions of the five targets employed in our second group of experiments and their estimates produced by all the considered algorithms are shown; note that, unlike FFT-based and MUSIC-based algorithms, the RASCAs achieve good accuracy even in the presence of closely spaced targets, like T 4 and T 5 .

XI. CONCLUSION
In this manuscript, six novel algorithms, dubbed range & angle serial cancellation algorithms (RASCAs), have been developed for the detection and the estimation of multiple targets in colocated MIMO radar systems. All these algorithms can be seen as instances of a general approach to target detection and estimation, and exploit new methods for the estimation of multiple overlapped real and complex tones. Moreover, four of them apply to colocated MIMO FMCW radar systems, whereas the other two to colocated MIMO SFCW radar systems. As evidenced by our computer simulations run on both synthetically generated data and measurements acquired through commercial devices, the devised algorithms are able to generate accurate 2D and 3D radar images in the presence of multiple closely spaced targets and outperform other algorithms based on the computation of multiple FFTs or on the MUSIC for DOA estimation.

APPENDIX COMPUTATIONAL COMPLEXITY
In this Appendix, the computational complexity, in terms of flops, is assessed for the RASCA-FC3 developed in Paragraph IV. The overall computational cost of this algorithm can be expressed as C FC3 = N VR C T 1 + N A K T 2 C T 2 + L b K T 2 C T 3 + C T sc , (191) where C T 1 is the contribution due to the first task of the RASCA-FC3, K T 2 (K T 3 ) represents the overall number of iterations carried out by the STDREC (STDAEC) algorithm, C T 2 (C T 3 ) is the contribution due to a single iteration of the STDREC (STDAEC) executed on a single VA (on the whole virtual array for a given frequency bin) and C T sc is the contribution due to the computation of the spatial coordinates of the overall image. The general criteria adopted in estimating the computational costs appearing in the RHS of (191) are illustrated in [38] and can be summarised as follows: T 1 -The cost C T 1 can be expressed as where: a) C x ZP = 4 N is the contribution due to the computation of the vectors {x (v) k,ZP ; k = 0, 1, 2} (see (43)-(45)); b) C X = 24N 0 log 2 N 0 is the contribution due to the computation of the vectors {X (v) k ; k = 0, 1, 2} (see (47)). T 2 -The computational cost of this task is mainly due to its main algorithm, i.e., to the STDREC algorithm. The cost C T 2 can be expressed as where: a) C CSFE = 4N CSFE I 2 is the cost originating from the CSFE 19 employed in STDREC-S1; b) C C X k = 18N 0 is the contribution due to the computation of the vectors (C ) (see (140)-(142)); c) C E = 8N 0 − 2 is the contribution due to the computation of the residual energy (see (54)).
T 3 -The cost C T 3 can be expressed as 19 Note that, in this case, the cost of the CSFE does not account for the evaluation of three DFTs, since these have been already evaluated in T1.
where: a) C CSFE V = 8N 0 log 2 (N 0 )+4N CSFE I 2 is the cost originating from the CSFE employed in STDAE-S1; b) C X (VF) = 6N VV N VH + 2N VV is the contribution due to the computation of the vertically folded spectrum X (VF) i [l] (76) in STDAE-S2; c) C CSFE H = 8N 0 log 2 (N 0 ) + 4N CSFE I 2 is the cost originating from the CSFE employed in STDAE-S3; d) C X OF = 6N VV N VH + 18N VV N VH N 0 is the contribution due to the computation of the overall folded spectrum {X m,OF [l]; m = 0, 1, 2} (see (78) and (83)); e) C CSFE OF = 4N CSFE I 2 is the cost due to the CSFE 20 in STDAE-S4; f) C E = 6N VV N VH is the contribution due to the computation of the residual energy in STDAEC-S3 (see (63)).
Finally, the cost C T sc = 5L is required to generate the overall point cloud.
Based on the results illustrated above, (191) can be rewritten as C FC3 = N VR (4N + 24N 0 log 2 (N 0 )) (195) where N T 2 N A K T 2 and N T 3 L b K T 2 .

ACKNOWLEDGMENT
The author would like to thank the anonymous Reviewers for their constructive comments, that really helped us to improve the quality of this manuscript.