Novel Methods for Approximate Maximum Likelihood Estimation of Multiple Superimposed Undamped Tones and Their Application to Radar Systems

In this manuscript, novel methods for the detection of multiple superimposed tones in noise and the estimation of their parameters are derived, and their application to colocated multiple-input multiple-output radar systems is investigated. These methods are based on a maximum likelihood approach and combine an innovative single tone estimator with a serial cancellation procedure. Our numerical results lead to the conclusion that these methods can achieve a substantially better accuracy-complexity tradeoff than various related techniques available in the technical literature. Moreover, they can be exploited to detect multiple closely spaced targets and estimate their spatial coordinates in radar systems.


I. INTRODUCTION
T HE problem of estimating the amplitude, phase and frequency of multiple (say, L) tones in additive white Gaussian noise (AWGN) has received significant attention for a number of years because of its relevance in various fields, including radar systems and wireless communications [1], [2]. It is well known that the maximum-likelihood (ML) approach to this problem leads to a complicated nonlinear optimization problem. Substantial simplifications can be made when the L tone frequencies are sufficiently well separated and the number N of available signal samples is large enough [3]- [5]. In fact, under these assumptions, each tone has a limited influence on the estimation of the others, so that approximate ML estimation can be achieved through a conceptually simple sequential procedure, that consists in iteratively executing two steps [4]. In the first step of this procedure, the parameters of the dominant tone (i.e., of the tone associated with the largest peak in the periodogram of the observed signal) are estimated in a ML fashion. In its second step, instead, the estimated tone is subtracted from the available signal samples and a new periodogram is computed for the resulting residual. These steps are repeated until all the detectable tones are estimated. Note that this simple sequential procedure can be considered as a specific instance of: a) the serial interference cancellation (SIC) concept [6], since, when a new tone is detected and its parameters estimated, all the other tones are considered as a form of interference; b) the so called CLEAN method [7], developed for high-resolution radio interferometry and later used in microwave imaging [8] and radar systems [9], [10]. Moreover, its technical relevance is motivated by the following relevant advantages [11]: 1. It turns a complicated multidimensional problem (whose dimensionality is usually unknown a priori) into a sequence of lower dimensional subproblems. Consequently, its overall complexity is proportional to that required to solve each of such subproblems and is usually much lower than that of: a) parametric estimation methods, like the MUSIC [12], the ESPRIT [13] and their computationally efficient versions (e.g., see [14] and [15]); b) non parametric spectral estimators, like the Capon method [16], the amplitude and phase estimation of a sinusoid (APES) [17], the iterative adaptive approach for amplitude and phase estimation (IAA-APES) [18] and their computationally efficient versions (e.g., see [19] and [20]).
2. It performs better than independently estimating the tones associated with the largest peaks of the original periodogram. In fact, it allows to identify peaks that are initially masked by the leakage due to nearby stronger tones.
3. It is able to estimate an unknown L in a simple fashion. In fact, this result can be achieved setting the initial value of this parameter to zero and applying a suitable test to establish whether, at each repetition of the first step, the largest peak detected in the periodogram of the last residual is significant [3] or whether, at each repetition of the second step, the energy of the new residual is large enough [21]. If one of these conditions is satisfied, the estimate of L is incremented by one and, then, the next step is carried out; otherwise, the estimation process is terminated. It is worth stressing that various estimation methods (e.g., the MUSIC and the ESPRIT estimators) require prior knowledge of L and that, in these cases, the use of some methods, like the generalized Akaike information criterion [22] or the minimum description length [23] (see also refs. [12] and [13]) is commonly proposed for the estimation of this parameter; however, the computational effort they require is not negligible.
The two-step procedure described above, despite its advantages, suffers from the following two shortcomings: 1. Any inaccuracy in the estimation of each single tone accomplished in its first step results in an imperfect cancellation of the tone itself and, consequently, in error accumulation; the intensity of this phenomenon increases with iterations, so affecting the estimation accuracy of the weakest tones.
2. The estimate of each tone is potentially biased, because of the presence of other tones [5]. Biases are influenced by the relative phase, frequency and amplitude of the superimposed tones and are expected to be more relevant in the first estimated frequencies, since these suffer from stronger interference from other tones. For this reason, the overall accuracy of this procedure depends on that of the employed single tone estimator and can be improved by adopting specific methods for mitigating the estimation bias. As far as the first issue is concerned, it is important to point out that optimal (i.e., ML) estimation of a single tone in AWGN is a computationally hard task. This is mainly due to the fact that the ML metric is an highly nonlinear function, that does not lend itself to easy maximisation (e.g., see [24]). In practice, the most accurate ML-based single tone estimators available in the technical literature achieve approximate maximisation of this metric through a two-step procedure; the first step consists in a coarse search of tone frequency, whereas the second one in a fine estimation generating an estimate of the so called frequency residual (i.e., of the difference between the real frequency and its coarse estimate). Coarse estimation is always based on the maximization of the periodogram of the observed signal, whereas fine estimation can be accomplished in an open loop fashion or through an iterative procedure. On the one hand, all the open loop estimators exploit spectral interpolation to infer the frequency residual from the analysis of the fast Fourier transform (FFT) coefficients at the maxima of the associated periodogram and at frequencies adjacent to it [11], [25]- [34]. Unluckily, unlike iterative estimators, the accuracy they achieve is frequency dependent and gets smaller when the signal frequency approaches the center of one of the FFT bins. On the other hand, the iterative estimation techniques available in the technical literature are based on various methods, namely on: a) standard numerical methods for locating the global maximum of a function (e.g., the secan method [35] or the Newton's method [36]); b) an iterative method for binary search, known as the dichotomous search of the periodogram peak [37]- [39]; c) interpolation methods amenable to iterative implementation [40]- [46]; d) the combination of the above mentioned dichotomous search with various interpolation techniques [47]; e) the computation of the first derivative of the spectrum [48].
The use of some of these algorithms in multiple tone estimators based on the above mentioned serial cancellation approach has been investigated in refs. [11], [21], [22], [49], [50] and [51]. More specifically, on the one hand, the periodogram-based (coarse) estimation method has been employed in the CLEAN algorithm [7], in the more CLEAN (MCLEAN) [21] and in the RELAX algorithm [22]. Note that, since a fine estimation step is missing in all these algorithms, achieving high accuracy requires the use of zero padding and of a large FFT order. On the other hand, the exploitation of more refined single tone estimators has been investigated in refs. [11], [49], [50] and [51]. In particular, the use of open-loop interpolation methods exploiting three or five adjacent spectral coefficients (including the one associated with the coarse frequency estimate) has been studied in ref. [11], whereas that of the iterative methods developed in refs. [40] and [52] has been analysed in refs. [49], [50] and [51], respectively.
As far as the second technical issue (i.e., estimation bias) is concerned, it is worth mentioning that the most straightforward methods for bias mitigation rely on the use of a) zero-padding for enhancing periodogram spectral resolution and b) window functions [5], [24], [53]- [56]; the price to be paid for these choices is an increase in the overall computational cost and in the variance of computed estimates, respectively. More refined methods are represented by interpolators with intrinsic leakage rejection [11] and nonlinear optimization methods. The last class of methods includes the expectation maximization (EM) algorithm [57], the space-alternating generalized EM (SAGE) algorithm [58], [59], the Newton's method [36], [60] and different optimization algorithms that employ cyclic cancellation procedures [21], [22], [51], [60]. In the last case, tone re-estimation is accomplished after removing the interference of both stronger and weaker tones as the iterations of the serial cancellation procedure evolve [21], [22] or after detecting and estimating the parameters of all tones [51]; the most refined version of the first method is described in ref. [22], where tone re-estimation is iterated after the estimation of each new tone, in order to generate excellent initial estimates for the next step (i.e., for the estimation of the next tone). Tone re-estimation reduces error accumulation and leads to convergence to the ML solution in the absence of noise if the frequency spacing of the detected tones is large enough; however, this result is achieved at the price of an increase of the overall computational cost and latency [21].
This manuscript aims at providing various new results about the estimation of multiple superimposed tones and their exploitation in colocated frequency modulated continuous wave (FMCW) and stepped frequency continuous wave (SFCW) multiple-input multiple-output (MIMO) radar systems operating at millimeter waves [61], [62]. Its contribution is threefold. First, a novel ML-based iterative estimator of a single real and complex tone is developed. The derivation of this estimator is based on: a) expressing the dependence of the ML metric on the frequency residual in an approximate polynomial form through standard approximations of trigonometric functions; b) exploiting the alternating minimization technique for the maximization of this metric (e.g., see [63, Par. IV-A]). Moreover, its most relevant feature is represented by the fact that it requires the evaluation of spectral coefficients that are not exploited by all the other related estimation methods available in the technical literature. Secondly, it is shown how serial cancellation in the frequency domain can be combined with our iterative estimator to detect multiple tones and estimate their parameters. Thirdly, the accuracy of our single and multiple estimators is assessed by extensive computer simulations; both synthetically generated data and the measurements acquired through two commercial MIMO systems in different scenarios are processed. Our results lead to the conclusion that our estimators outperform all the other related estimators in terms of probability of convergence and accuracy in the presence of arbitrary frequency residuals.
The remaining part of this manuscript is organised as follows. In Section II, signal models are defined and their relevance in colocated MIMO radar systems is discussed. Section III is devoted to the derivation of our single tone and multiple tone estimation algorithms, to the assessment of their computational complexity, and to analysis of their similarities and differences with related estimators available in the technical literature. In Section IV, the performance of our estimation algorithms is assessed and compared with that achieved by other estimators. Finally, some conclusions are offered in Section V.

II. SIGNAL AND SYSTEM MODELS
In this manuscript, we focus on the problem of estimating all the parameters of the complex sequence and its real counterpart with n = 0, 1, ..., N − 1; here, A l and F l ∈ [0, 1) denote the complex amplitude and the normalised frequency, respectively, of the l-th complex tone appearing in the right-hand side (RHS) of eq. (1), a l |A l |, ψ l ∠(A l ), w c,n is the n-th sample of an additive white Gaussian noise (AWGN) sequence (whose elements have zero mean and variance 2σ 2 ), w r,n {w c,n }, N is the overall number of samples, and {x} and arg(x) denote the real part and the phase, respectively, of the complex quantity x. It is useful to point out that the signal model (2) can be rewritten as where represents the complex amplitude of the real tone appearing in the RHS of eq. (2). The signal models (1) and (2) appear in a number of problems concerning radar systems, wireless communications and biomedical applications. In the remaining part of this paragraph, we analyse their meaning in the context of colocated and bistatic multiple-input multiple output (MIMO) radar systems operating in time division multiplexing (TDM) [1] and at millimetre waves. These systems have the following important characteristics: a) They are equipped with a transmit (TX) antenna array and a receive (RX) antenna array, that consist of N T elements and N R elements, respectively; these allow to radiate orthogonal waveforms from different antennas and to receive distinct replicas of the electromagnetic echoes generated by multiple targets. b) Their TX and RX arrays are made of distinct antenna elements, placed at different positions. Moreover, TX antennas are close to the RX ones and, in particular, are usually placed on the same shield. c) Orthogonality of the transmitted waveforms is achieved by radiating them through distinct TX antennas over disjoint time intervals.
The architecture of the radar system considered in this manuscript is illustrated in Fig. 1. In the remaining part of this section, two different models are described for the RF signal generated by the voltage controlled oscillator (VCO) of its transmitter and radiated by its TX array after power amplification. In the first case, corresponding to a frequency modulated continuous wave (FMCW) radar system, the VCO is fed by a periodic ramp generator; this produces a chirp FM signal, whose instantaneous frequency evolves periodically, as illustrated in Fig. 2. 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 [1], whereas the parameters f 0 and B are the start frequency and the bandwidth, respectively, of the transmitted signal. If we assume that the radar system exploits all the available TX diversity (i.e., all its TX antennas) and that a time slot of T 0 s is assigned to each TX antenna, transmission from the whole TX array is accomplished over an interval lasting T F N T T 0 s; this interval represents the duration of a single  Fig. 2: Representation of the instantaneous frequency of the RF signal generated by the VCO in a FMCW radar system. transmission frame. 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, s(t) exp(jθ(t)), and µ B/T is the chirp rate (i.e., the steepness of the generated frequency chirp). Let r (q) RF (t) denote the signal available at the output of the q-th RX antenna in the same time interval, with q = 0, 1, ..., N R − 1 (see Fig. 1). This signal feeds a low noise amplifier (LNA), whose output undergoes downconversion, filtering and analog-to-digital conversion at a frequency f s = 1/T s ; here, T s denotes the sampling period of the employed analog-to-digital converter (ADC). If the radiated signal s RF (t) (5) is reflected by L static point targets, the useful component of r RF (t) consists of the superposition of L echoes, each originating from a distinct target. In this case, if the propagation environment is static or undergoes slow variations, a simple mathematical model can be developed to represent the sequence of samples generated by the ADC in a single chirp interval. In deriving that model, the couple of the involved physical TX and RX antennas (namely, the p-th TX antenna and the q-th RX antenna) of the considered bistatic radar is usually replaced by a single virtual antenna of an equivalent monostatic radar. In particular, if it is assumed that all the TX and RX antennas belong to the same planar shield and a reference system lying on the physical antenna array is defined, the abscissa x v and the ordinate y v of the v-th virtual antenna element associated with the p-th TX antenna and the q-th RX antenna can be computed as (e.g., see [64,Par. 4 and respectively, with v = 0, 1, ..., N V R − 1; here, (x p , y p ) ((x q , y q )) are the coordinates of the TX (RX antenna) and N V R N T · N R represents the overall size of the resulting virtual array. Based on these assumptions, the n-th received signal sample acquired through the v-th virtual antenna element (with v = 0, 1, ..., N V R − 1) can be expressed as (e.g., see [65,Par. 4.6,eq. (4.27)]) with n = 0, 1, ..., N − 1; here, N is the overall number of samples acquired over a chirp period, a l is the amplitude 1 of the l-th component of the useful signal, F is the normalized version of the frequency f characterizing the l-th target detected on the v-th virtual RX antenna, 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, respectively, and w (v) r,n , the n-th sample of the AWGN sequence affecting the received signal, is a Gaussian random variable having zero mean and variance σ 2 (assumed to be independent of v). It is important to point out that: a) the real signal model (10) is identical to that expressed by eq. (2) and can be adopted in all the FMCW radar systems that extract 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 [62, Par 2.1 eq. (2.2)]). In the last case, the complex model equivalent to that expressed by eq. (1), must be adopted in place of its real counterpart (10) for any n; here, for any v and l.
The second case we consider in the generation of the radiated waveform corresponds to a stepped frequency continuous wave (SFCW) radar system. Its name is motivated by the fact that the VCO of its transmitter is fed by a staircase generator. For this reason, the instantaneous frequency of the resulting RF signal takes on N distinct and uniformly spaced values in an interval lasting T s for each TX antenna; the n-th value of the instantaneous frequency is with n = 0, 1, ..., N − 1; here, f 0 is the minimum radiated frequency, ∆f is the frequency step size and N is the overall number of transmitted frequencies. It can be shown that, under the same assumptions made in the derivation of eq. (15), the measurement acquired through the v-th virtual element at the n-th frequency can be expressed as with v = 0, 1, ..., N V R − 1; here, the amplitude ψ  are still expressed by eqs. (14) and (16), respectively, is the normalised delay characterizing the l-th target and observed on the v-th virtual antenna, and 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 (15). Then, in both the considered FMCW and SFCW radar systems, the signal observed on each virtual channel can be represented 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 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 eqs. (10)-(12) ((15)- (16) and (18)- (19)), 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 [66]. In fact, if, in a FMCW radar system, the l-th tone is found at the frequencyf l , the presence of a target at the range (see eqs. (12) and (13) is detected. Similarly, in a SFCW radar system, the normalised delayF (v) l estimated on the the v-th virtual channel is associated with a target whose range is (see eqs. (13) and (19) 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 V R phases {ψ (v) l ; v = 0, 1, ..., N V R − 1} observed over the available virtual antennas. In fact, since (see eqs. (13) and (14)) where is the wavelength associated with the frequency f 0 , the sequence {ψ 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 FFT processing on the estimated phases taken across multiple elements of the virtual array in a single frame interval [67], [68]. From the considerations illustrated above, the following conclusions can be easily inferred: a) on the one hand, achieving precise estimation of the range of a given target requires the availability of an accurate estimate of the normalised frequency (or delay) of the real or complex tone associated with the target itself over at least one RX antenna; b) on the other hand, the quality of the estimate of the direction of arrival (DOA) characterizing the radar echo generated by a given target depends on the accuracy of the phase estimated over multiple virtual channels. The last consideration motivates the importance of accurately estimating this parameter over multiple antennas. It is worth stressing that limited attention is often paid to this technical issue in most of the technical literature dealing with the estimation of the parameters of a single tone or multiple superimposed tones; in fact, a lot of emphasis is usually put on the accuracy of frequency estimation, but the quality of phase estimates is neglected (e.g., see [11], [21], [40], [49]- [51]). Finally, in analysing the suitability of multiple tone estimators to colocated MIMO radar systems operating at millimetre waves, the following additional 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 of resulting radar images is a cloud of point targets whose mutual spacing can be very small. For this reason, the accuracy of these images depends on the frequency resolution (delay resolution) achieved on each virtual antenna in a FMCW (SFCW) radar system. In fact, this makes the receiver able to separate point targets characterized by similar ranges.
2) Information about the RCS of each point targets can be exploited in the classification of extended targets; for this reason, the availability of an accurate estimate of the amplitude of each radar echo can be very useful in a number of applications (e.g., in SAR imaging [17], [22]).
3) Distinct radar echoes can be characterized by substantially different signal-to-noise ratios (SNRs). 4) 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 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. APPROXIMATE MAXIMUM LIKELIHOOD ESTIMATION OF SINGLE AND MULTIPLE TONES
In this section, we first derive a new method for estimating the parameters of a single real or complex tone. Then, we show how this method can be exploited to detect multiple superimposed tones and estimate their parameters through a deterministic procedure based on successive cancellation. Finally, we analyse the computational complexity of the developed estimation methods and compare them with some related techniques available in the technical literature.

A. Estimation of a single tone
Let us focus first on the problem of estimating the parameters (namely, the frequency and complex amplitude) of a single tone contained in the real sequence {x r,n ; n = 0, 1, ..., N − 1}, whose n-th sample is expressed by eq. (10) with L = 1, i.e. as x r,n = a cos(2πnF + ψ) + w r,n (27) or, equivalently, as (see eq. (3)) where (see eq. (4)) and F are the complex amplitude and the normalised frequency, respectively, of the tone itself. It is well known that the ML estimates F ML and C ML of the parameters F and C, respectively, represent the solution of the least square problem (e.g., see whereF andC represent trial values of F and C, respectively, is the mean square error (MSE) evaluated over the whole observation interval, is the square error between the noisy sample x r,n (28) and its useful component evaluated under the assumption that F =F and C =C. The cost function ε(F ,C) (31) can be easily reformulated in a different way as follows. To begin, substituting the RHS of eq. (33) in that of eq. (32) produces, after some manipulation, and {x} ( {x}) denotes the real (imaginary) part of the complex variable x. Then, substituting the RHS of eq. (34) in that of eq. (31) yields where is, up to the scale factor 1/N , the Fourier transform of the sequence {x r,n }, Based on eq. (36), it is not difficult to show that the optimization problem expressed by eq. (30) does not admit a closed form solution because of the nonlinear dependence of the function ε(F ,C) on its variableF . However, an approximate solution to this problem can be derived by a) Exploiting an iterative method, known as alternating minimization (AM) [63]. This allows us to transform the twodimensional (2D) optimization problem expressed by eq. (30) into a couple of interconnected one-dimensional (1D) problems, one involving the variableF only (conditioned on the knowledge ofC), the other one involving the variableC only (conditioned on the knowledge ofF ). b) Expressing the dependence of the function ε(F ,C) on the variableF through the couple (F c ,δ) such that where F c is a given coarse estimate of F ,δ is a real variable called residual, is the normalized fundamental frequency associated with the N 0 -th order discrete Fourier transform (DFT) of the zero padded version of the vector collecting all the elements of the sequence {x r,n }, M is a positive integer (dubbed oversampling factor), 0 D is a D−dimensional (column) null vector and N 0 M · N . c) Expressing the dependence of the function ε(F ,C) (36) on the variableδ through its powers {δ l ; 0 ≤ l ≤ 3}; this result is achieved by approximating various trigonometric functions appearing in the expression of ε(F ,C) with their Taylor expansions truncated to a proper order.
Let us show now how these principles can be put into practice. First of all, the exploitation of the above mentioned AM approach requires solving the following two sub-problems: P1) minimizing the cost function ε(F ,C) (36) with respect toC, givenF =F ; P2) minimizing the same function with respect toF , givenC =Ĉ. Sub-problem P1 can be easily solved thanks to the polynomial dependence of the cost function ε(F ,C) on the variablesC R andC I . In fact, the minimum of the function ε(F ,C) with respect to the last variables is the solution of the equations andC that result from computing the partial derivative of the RHS of eq. (36) with respect toC R andC I , respectively. Solving the linear system of equations (45)- (46) in the unknownsC R andC I produces the optimal valueŝ ofC R andC I , respectively. Putting together the last two formulas yieldŝ Therefore, givenF =F , the optimal valueĈ of the variableC can be computed exactly through the last equation; this requires the evaluation ofX(F ) and g(F ) (see eqs. (38) and (39), respectively). Note that, on the one hand, g(F ) can be easily evaluated through its exact expression which is easily derived from its definition (39). On the other hand,X(F ) can be computed exactly through its expression (38) or, in an approximate fashion, through a computationally efficient procedure based on the fact that the vector collects N 0 uniformly spaced samples of the functionX(F ) (38). In fact, the k-th element of the vector X 0 (42) can be expressed asX where Let us take into consideration now sub-problem P2. Unluckily, this sub-problem, unlike the previous one, does not admit a closed form solution. For this reason, an approximate solution is developed below. Such a solution is based on representing the normalized frequency F in the same form asF (see eq. (40)), i.e. as and on a novel method for estimating the real residual δ, i.e. for accomplishing the fine estimation of F . This method is derived as follows. Representing the trial normalized frequencyF according to eq. (40) allows us to express the variableφ n (35) as Then, substituting the RHS of eq. (56) in that of eq. (34) (withC =Ĉ) yields Based on the last equation and some trigonometric approximations (see Appendix A), the approximate expression can be derived for the function ε(F ,Ĉ) (36); here, for any ρ and k = 1, 2, and x k,n n k · x r,n with n = 0, 1, ..., N − 1. It is important to point out that: a) if ρ is an integer, the quantityX k,ρ (63) (with k = 1 and 2) represents the ρ-th element of the vector generated by the N 0 -th order DFT of the zero padded version of the vector b) if ρ is not an integer, the quantityX k,ρ can be evaluated exactly on the basis of eq. (63) or, in an approximate fashion, by interpolating a subset of adjacent elements of the vector X k (65); c) the evaluation of the vectors {X k ; k = 1, 2} requires two additional FFTs.
Since the function ε SFE (∆,Ĉ) (60) is a polynomial of degree three in the variable∆, an estimate∆ of ∆ and, consequently, an estimate (see eq. (57) of δ, can be obtained by computing the derivative of this function with respect to∆, setting it to zero and solving the resulting quadratic equation in the variable∆; here, and Note that only one of the two solutions of eq. (69), namelŷ has to be employed. A simpler estimate of ∆ is obtained neglecting the contribution of the first term in the left-hand side of eq. (69), i.e. setting a(ρ) = 0. This leads to a first-degree equation, whose solution iŝ Given an estimate∆ of ∆ (and, consequently, and estimateδ of δ; see eq. (68)), the fine estimatê of F can be evaluated on the basis of eq. (55). The mathematical results derived above allow us to derive a novel estimation algorithm, called single frequency estimator (SFE), for iteratively estimating the normalised frequency F and the complex amplitude C. This algorithm is initialised by 1) Evaluating: a) the vector X 0 (42); b) the initial coarse estimateF where the integerα is computed by means of the well known periodogram method (e.g., see [35,Sec. IV] c) the quantity (see eq. (61))ρ 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 is 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).
At the end of the last (i.e., of the N SFE -th) iteration, the fine estimatesF =F (NSFE) andĈ =Ĉ (NSFE) of F and C, respectively, become available.
The SFE is summarized in Algorithm 1. It is important to point out that: a) The estimateδ (i) of the residual δ computed by the SFE in its i-th iteration is expected to become smaller as i increases, sinceF (i) should progressively approach F if our algorithm converges.
b) The coefficients {K p (2ρ); p = 1, 2, 3} can be computed exactly on the basis of eq. (62). However, since the definition (62) can be put in the equivalent form where the identities and holding for any q ∈ C, can be exploited for the efficient computation of K p (2ρ) for any p and ρ.
c) The quantities {X k,ρ (i−1) ; k = 1, 2} required in the first step of the i-th iteration can be computed exactly on the basis of eq. (63). However, they can be also evaluated, in an approximate fashion, by interpolating I adjacent elements of the N 0 -dimensional vectors X k (65), where I denotes the selected interpolation order.
d) The estimate∆ (i) evaluated according to eq. (74) is expected to be less accurate than that computed on the basis of eq. (73). However, our numerical results have evidenced that both solutions achieve similar accuracy. Despite, the use of eq. (73) is assumed in Algorithm 1 for generality.
e) The SFE can be employed even if the single tone appearing in the RHS of eq. (27) is replaced by the superposition of L distinct tones (see eq. (3)). In this case, the strongest (i.e., the dominant) tone is detected through the periodogram method (see eq. (77)) and the parameters of this tone are estimated in the presence of both Gaussian noise and the interference due to the remaining tones. Therefore, the estimation accuracy of the SFE is affected by both the amplitudes and the frequencies of the other (L − 1) tones.
f) A stopping criterion, based on the trend of the sequence {∆ (i) ; i = 1, 2, ...}, can be easily formulated for the SFE. For instance, the execution of its two steps can be stopped if, at the end of the i-th iteration, the condition is satisfied; here, ε ∆ represents a proper threshold. All the results developed above refer to the real sequence {x r,n }, whose n-th element is expressed by eq. (27). However, a similar estimation method (dubbed complex SFE, CSFE) can be developed for the complex counterpart, i.e. for a complex sequence {x c,n } such that (see eq. (15) with L = 1) with n = 0, 1, ..., N − 1. In this case, the ML estimation of the parameters F and A can be formulated in a similar way as eq. (30), the only differences being represented by the fact that: a) the parameter C is replaced by A; b) the term ε n (F ,C) appearing in the RHS of eq. (31) is replaced by where represents the useful component of x c,n (89) under the assumption that F =F and A =Ã. Substituting the RHS of the last equation in that of eq. (90) and, then, the resulting expression in the RHS of the MSE (see eq. (1)) yields, after some manipulation, where The procedure adopted to minimise the last function with respect toF andÃ is similar to that illustrated above for the SFE. For this reason, in the following we limit to show essential mathematical results only. First all, it is easy to show that, giveñ F =F , the function ε(F ,Ã) (93) is minimised with respect toÃ selecting 2 whereX(F ) can be computed exactly through its expression (see eq. (38)) replacing the real sequence {x r,n } with the complex sequence {x c,n }) or, in an approximate fashion, through an interpolation technique following the same mathematical approach as that illustrated for the SFE. Note also that the mathematical result expressed by eq. (95) is exact. On the contrary, giveñ A =Â, a closed form expression for the value ofF minimising the function ε(F ,Â) cannot be derived because of the nonlinear dependence of this function onF . However, following the same mathematical approach as that illustrated for the SFE, the approximate expression can be obtained for the cost function ε(F ,Ã) (93) (see Appendix A for further details); the quantities appearing in the last equation have the same meaning as the ones defined for eq. (60). It worth stressing that: a) The considerations expressed about the evaluation of the quantities {X k,ρ } in our derivation of the SFE apply to the CFSE too. However, in this case, three additional spectral coefficients, namely {X k,ρ ; k = 1, 2, 3}, need to be computed.
b) The initial estimateρ (0) =α (see eq. (78)) of ρ is evaluated in a similar way as the SFE (see eq. (77)), i.e., aŝ α = arg max α∈{0,1,...,N0−1} The minimization of the function ε CSFE (F ,Ã) (96) with respect to∆ is achieved by taking its partial derivative with respect to this variable and setting it to zero; this results again in the quadratic equation (69), whose coefficients are and Following the approach adopted in the development of the SFE and exploiting the mathematical results expressed by eq. (95) and by eqs. (98)-(100) allow us to easily derive the CSFE; this algorithm is summarized in Algorithm 2. It is worth noting that, if the number of iterations is sufficient, the SFE and the CSFE algorithms are asymptotically unbiased.

B. Estimation of multiple tones
Let us analyse now in detail how the techniques derived in the previous paragraph can be exploited to estimate the multiple tones that form the useful component of the real input sequence {x r,n } (the complex sequence {x c,n }), when its n-th sample is expressed by eq. (2) (eq. (1)) with L > 1. The two recursive methods we develop to achieve this target are based on the following basic principles: 1) Tones are sequentially detected and estimated.
2) The detection of a new tone and the estimation of its parameters are based on the procedure developed for the SFE and the CFSE in the previous paragraph; in addition, a cancellation algorithm is incorporated in both these methods to remove the contribution of previously detected tones from all the spectral information (namely, the spectrumX(F ) (38), the vector X 0 (42) and the coefficients {X k,ρ } (63)), that are processed to detect and estimate the new tone. 3) After detecting a new tone and estimating its parameters, a re-estimation technique is executed to improve the accuracy of both this tone and the previously estimated tones; the proposed technique is inspired by the related methods described in refs. [11], [21] and [22]. 4) A proper criterion is adopted to stop recursions. This allows to estimate the (unknown) number of targets, that is the value of the parameter L. The first recursive method relying on these principles is called single frequency estimation and cancellation (SFEC) or complex single frequency estimation and cancellation (CSFEC), depending on the fact that its input sequence is real or complex, respectively. In the following, for space limitations, we focus on the CSFEC method only; however, readers should keep in mind that the differences between this method and the SFEC are similar to those illustrated in our description of the single tone estimation algorithms on which they are based.
The CSFEC algorithm is initialised by: 1) Executing the CSFE, fed by the complex sequence {x c,n }, to generate, through N CSFE iterations, the initial estimateŝ F 0 [0] andÂ 0 [0] of the normalized frequency and the complex amplitude, respectively, of the first detected tone.
2) Setting the recursion index r to 1. Then, a recursive procedure is started. The r-th recursion is fed by the vectorŝ collecting the frequencies and the associated complex amplitudes characterizing the r tones detected and estimated in the previous (r − 1) recursions, and generates the new vectorŝ after: a) estimating the frequencyF r [r] and the associated complex amplitudeÂ r [r] of a new (i.e., of the r-th) tone (if any); b) refining the estimates of the r tones available at the beginning of the considered recursion. The procedure employed for accomplishing all this consists of the three steps described below (the p-th step is denoted CSFEC-Sp). CSFEC-S1 (CFSE initialisation with cancellation) -In this step, the following quantities are evaluated (see the initialisation part of Algorithm 2): a) The residual spectrum where X 0 is the N 0 -th order DFT of the zero padded version x 0,ZP of the vector x 0 collecting all the elements of the sequence {x c,n } (see eqs. (43)-(44)) and the N 0 -dimensional vector represents the contribution given by all the estimated tones to X 0 (in particular, ) is the contribution provided by l-th tone to the vector X 0 (the expression of this vector can be found in Appendix B). If the overall energy satisfies the inequality where T CSF EC is a proper threshold, the algorithm stops and the estimateL = r of L is generated. b) The integer (see eq. (97))α [r] = arg max α∈{0,1,...,N0−1}  Our simulation results have evidenced that, unluckily, the estimates generated by the SFEC and CSFEC algorithms are biased. However, this bias can be removed by running an additional (i.e., a fourth) step after that these algorithms have been executed. In this final step, the estimation algorithm developed by Ye and Aboutanios in ref. [49], [50]  It is worth pointing out that: a) The oversampling factor M adopted in the computation of the vectors {X (l) k } and the stopping criterion employed by the SFE (or CSFE) need to be carefully adjusted in order to achieve a good accuracy in the estimation of the parameters of each new tone. b) Poor estimation of the normalised frequency F l and/or the complex amplitude A l (C l ) may lead to significant error accumulation if CSFEC-S3 (the corresponding step in the SFEC) is removed; readers should always keep in mind that a fundamental role in accurate cancellation is played by the accuracy of the estimated frequency residual.
c) The threshold T CSFEC needs to be properly adjusted in order to ensure that the probability thatL equals to L is close to unity. On the one hand, a large value of T CSFEC may lead to miss weaker tones; on the other hand, a small value of this parameter may lead to the identification of nonexistent tones and, consequently, of false targets in a radar system.

C. Estimation of single and multiple delays
All the estimation techniques developed in the previous two paragraphs refer to the real (complex) sequence {x r,n } ({x c,n }), whose n-th element is expressed by eq. (10) (eq. (15)). As shown in Section II, these models are suitable to represent the sampled downconverted signal made available by the receiver of a colocated FMCW MIMO radar. However, similar estimation methods can be also developed for a colocated SFCW MIMO radar system, since the complex signal model (18) illustrated for this system is very similar to the model (15) referring to its complex FMCW counterpart, the only differences being represented by the fact that: a) the physical meaning of the parameter F (v) l , since this represents a normalized delay in eq. (18) and a normalized frequency in eq. (15); b) the sign of the argument of all the complex exponentials appearing in the RHS of eq. (18), since this is the opposite of that of the corresponding functions contained in the RHS of eq. (15). The implications of this sign reversal can be summarised as follows. In the derivation of the single tone estimator, the functionX(F ) (38) is replaced byX where x c,n is the n-th sample of the received signal sequence (see eq. (18)). Consequently, the vectors X 0 (42) and X k (65) are generated by taking the N 0 order inverse DFT (IDFT) of the zero-padded vectors x 0,ZP (43) and x k,ZP (66), respectively; note that the n-th element of the N 0 -th dimensional vector x 0,ZP (x k,ZP ) is equal to x c,n (n k · x c,n ; see eq. (64)) for n = 0, 1, ..., N − 1 and is equal to zero for n ≥ N . An algorithm for estimating a single delay can be easily developed by following the same line of reasoning as the CSFE; this is called complex single delay estimator (CSDE). More specifically, the CSDE computes an estimateδ of the delay residual δ by solving eq. (69) in the variable∆ (57); in this case, the coefficients a(ρ) and c(ρ) are still expressed by eqs. (98) and (100), respectively, whereas whereÂ is the estimate of the complex gain A computed on the basis of eq. (95). Moreover, in the evaluation of the coefficients {a(ρ), b(ρ), c(ρ)} of eq. (69), the quantityX k,ρ (i−1) (see eq. (63)) is replaced bȳ whereρ (i−1) is still defined by eq. (80). These considerations must be kept into account in the development of the counterpart of the CSFEC algorithm; this is called complex single delay estimation and cancellation (CSDEC). Also in this case, the estimation algorithm developed by Ye and Aboutanios in ref. [49], [50] is interconnected to our CSDEC algorithm; the resulting hybrid version is dubbed hybrid CSDEC (HCSDEC).

D. Comparison with other estimation methods
The SFEC and CSFEC techniques developed in the Paragraphs III-A-III-B are conceptually related to the multiple tone estimators developed by Gough [21], Li and Stoica [22], Macleod [11], Ye and Aboutanios [49], [50], and Serbes and Qaraqe [51] (these algorithms are denoted CLEAN, RELAX, Alg-M, Alg-YA and CFH, 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. Despite their similar structure, they exhibit various differences, that concern the two specific issues listed below.
Single frequency estimator -The main difference is represented by the algorithm they employ in the estimation of a single tone. In fact, on the one hand, the CLEAN and RELAX algorithms rely on the coarse estimate generated by the periodogram method for each detected tone and, eventually, exploit zero-padding to improve spectral resolution. On the other hand, the Alg-M, the Alg-YA and the CFH algorithm compute the frequency residual by means of open-loop interpolation or iterative methods; the last methods refine the estimate of the frequency residual through multiple iterations. All the single tone estimators employed in these algorithms differ from the ones used in the SFEC and CSFEC, and their hybrid versions.
Use of a re-estimation procedure -In the CLEAN and RELAX algorithms and in Alg-M, once a new tone has been estimated, all the previously computed tones are re-estimated by subtracting the contribution of all the other tones; tone cancellation is accomplished in the time domain in the CLEAN and RELAX algorithms, whereas is carried out in the frequency domain in Alg-M. The last approach is also adopted in our estimation algorithms. Finally, the CFH algorithm and the Alg-YA accomplish re-estimation after computing a coarse estimate of all the parameters of the detected tones. However, the former algorithm executes this only once, whereas the latter one repeats it a given number of times.

E. Computational complexity
The complexity of the estimation algorithms developed in Paragraphs III-A-III-C has been carefully assessed in terms of number of floating operations (flops) to be executed in the detection of L targets. The general criteria adopted in estimating the overall computational cost of the SFE and that of the CSFE algorithms 4 are summarised in Appendix C, where a detailed analysis of the contributions due to the different tasks accomplished by each of them is also provided. Our analysis leads to the conclusion that these are approximately of order O(M SFE ) and O(M CSFE ), respectively, with  and if no re-estimation is accomplished (see CSFEC-S3 in the description of the CSFEC algorithms) and the algorithms stop after detecting the last tone. Note that the first term appearing in the RHS of both eqs. (122)-(123) accounts for the initialization (and, in particular, for the computation of the vectors X 0 (42) and {X k } (65)), whereas the second term for the fact that, in the SFEC (CSFEC), the SFE (CSFE) is executed L times. It is worth noting that the costs due to the evaluation of the estimated tones detected after the first one and to their frequency domain cancellation do not play an important role in this case. However, if re-estimation is accomplished, the parameter L appearing in the RHS of (122)-(123) is replaced by L 2 , since this task involves all the estimated tones. Despite this, the increase in the overall computational cost of the CSFEC (SFEC) with respect to the CSFE (SFE) is limited since, as evidenced by our simulation results, the use of re-estimation allows these algorithms to achieve convergence with a smaller value of the parameter N SFE (N CSFE ). Finally, it is important to compare the computational cost of the SFEC and CSFEC algorithms with that of the CLEAN, RELAX, Alg-M, Alg-YA and CFH techniques considered in the previous paragraph. Their order of complexity is listed in Table I, from which the following considerations can be easily inferred: 1) The CLEAN and the RELAX algorithm are characterized by the same order, espressed by the complexity of a zero-padded FFT multiplied by L 2 ; the last factor is due to the fact that tone re-estimation is employed in both algorithms.
2) The Alg-M is characterized by the lowest computational cost; in fact, since it performs the cancellation of the detected tones directly in the frequency domain, tone estimation does not require the computation of additional FFTs. Moreover, since tone re-estimation is employed, the cost for the estimation of a single tone (i.e., the parameter K M ) is multiplied by L 2 .
3) The order of complexity of the Alg-YA and that of the CFH algorithm depend on the fact that a FFT is computed for each tone; however, an additional term equal to L Q N is included in the order of Alg-YA since this estimates all the tones Q times.
4) The order of the cost of the SFEC (CSFEC) algorithm is similar to that of the Alg-M; however, in this case, three (four) FFTs are computed in the initialization phase and the cost for the estimation of each tone is multiplied by N SFE (N CSFE ) since tone refinement is performed within each recursion.

IV. NUMERICAL RESULTS
In this section, we compare, in terms of accuracy and computation time (CT), our single frequency estimator with the A&M [40], QSE and HAQSE algorithms [52], and our multiple tone estimator with the CFH algorithm [51] and the Alg-YA [49], [50]. As far as the A&M is concerned, two versions of it are considered; such versions are denoted A&M#1 and A&M#2 in the following and correspond to the Alg-1 and Alg-2, respectively, described in ref. [40]. Our numerical results are divided in two parts. The first part refers to the results generated on the basis of synthetically generated measurements, whereas the second one to the results generated on the basis of the measurements acquired through two commercial colocated MIMO radars. It is assumed that, in all the considered scenarios, a complex measurement sequence is observed; for this reason, the complex version of our algorithms is always taken into consideration.
A. Numerical results based on synthetically generated measurements Four different scenarios have been considered. The main features of these scenarios can be summarised as follows: Scenario #1 (S1) -This is characterized by L = 1, i.e. by a single tone, having amplitude equal to one and whose normalised frequency (phase) is uniformly distributed over the interval [F min , F min + F r ] ([0, 2π]), where N = 512 represents the overall number of samples, F min = 8.3/N and F r = 20/N . Moreover, the sampling frequency f s = N Hz has been selected in this scenario and in the remaining ones.
Scenario #2 (S2) -This is characterized by L = 2, i.e. by a couple of tones, whose amplitudes are equal to one and whose phases are uniformly distributed over the interval [0, 2π]. Moreover, their normalized frequencies are F 0 = 8.3/N and F 1 = 9.8/N , where N = 512 represents the overall number of samples.
Scenario #3 (S3) -This is characterized by L ∈ {2, 3,..., 10}, i.e. by a varying number of tones, whose amplitudes are equal to one and whose phases are uniformly distributed over the interval [0, 2π]. Their normalized frequencies, instead, are given by with k = 0, 1, ..., L − 1 and F 0 = 8.3/N . Scenario #4 (S4) -This is characterized by L = 10 tones, whose amplitudes decrease with their frequency and whose phases are uniformly distributed over the interval [0, 2π]. In particular, the amplitude of the k-th tone is given by with ∆ a = 7.5/(L − 1), whereas its normalized frequency is with F 0 = 8.3/N and k = 0, 1, ..., 9; here, represents the normalized frequency spacing between adjacent tones and is controlled through the non negative integer parameter m (m = 0, 1, ..., 5 is assumed in the following). In all the considered scenarios, the following performance indices have been assessed for each of the considered algorithms: a) the root mean square error (RMSE) for the normalized frequency estimates (RMSE f ) and for the phase estimates (RMSE p ); b) the failure probability (P f ), i.e. the probability that the considered algorithm does not converge. It is important to point out that: 1) in our simulation, a failure event is detected when the absolute value of the normalised frequency error is exceeds the threshold i.e., the size of half a frequency bin in the absence of oversampling; 2) the simulation runs in which an algorithm fails are ignored in the evaluation of its RMSEs; 3) each performance index is computed running 20000 Monte-Carlo simulations. Some numerical results referring to S1 are illustrated in Figs. 3, 4 and 5 (in all the following figures, simulation results are indicated by markers, whereas continuous lines are drawn to fit them, so facilitating the interpretation of the available data). More specifically, the dependence of RMSE f and RMSE p on the signal-to-noise ratio (SNR) is represented in Figs. 3 and 4, respectively; in both figures, the SNR range [−15, 20] dB is considered and the Cramer-Rao lower bound (CLRB) is also shown for comparison [70]. Five single tone estimators, namely the CSFE, A&M#1, A&M#2, QSE and HAQSE have been considered; moreover, the following parameters have been selected for them: a) overall number of iterations N CSFE = 20, interpolation order 5 I = 7, oversampling factor M = 4 for the CSFE; b) overall number of iterations N A&M = 2 for the A&M algorithms; c) N QSE = 3 for the QSE algorithm; d) N HAQSE = 2 for the HAQSE algorithm; e) q = N −1/3 for the QSE and HAQSE algorithms. In Fig. 5, instead, the dependence of P f on the SNR is shown for all these algorithms. From these results it is easily inferred that: 1) The CSFE technique appreciably outperforms all the other estimation techniques for any value of SNR both in frequency and phase estimation. 2) The CSFE accuracy approaches the CRLB both in frequency and in phase estimation; however, the RMSE of the frequency estimate slowly departs from the CRLB for large SNRs because of its imperfect convergence due to the limited number of iterations.
3) The CSFE technique is more robust than all the other estimators, since its failure probability is substantially lower than that of all the other algorithms for SNR > −10 dB.     Some numerical results referring to S2 are illustrated in Figs. 6, 7 and 8. More specifically, the dependence of RMSE f and RMSE p on the average SNR 6 is represented in Figs. 6 and 7, respectively; once again, the SNR range [−15, 20] dB has been considered. The accuracy of five multiple tone estimators, namely the CSFEC, HCSFEC and CFH algorithms, and the Alg-YA have been assessed; moreover, the following parameters have been selected for these algorithms in S2, S3 and S4: a) overall number of iterations N CSFE = 5, number of re-estimations N RES = 1, interpolation order I = 7, oversampling factor M = 4 for the CSFEC; b) number of Alg-YA iterations N Y A = 3 for the HCSFEC; c) the same parameters as the HAQSE in S1 for the CFH algorithm; d) overall number of iterations Q = 30 for Alg-YA. The dependence of P f on the SNR is illustrated in Fig. 8 for the same algorithms. These results show that: 1) The CFH algorithm is outperformed by all the other estimation techniques for any value of SNR both in frequency and phase estimation. 2) The CSFEC accuracy approaches the CRLB both in frequency and in phase estimation for SNR < −10 dB; for higher values of SNR, the error introduced by its frequency bias becomes visible.
3) The HCSFEC and the Alg-YA accuracy approach the CRLB both in frequency and in phase estimation for any value of SNR. 4) The CSFEC and the HCSFEC techniques are more robust than all the other estimators, since their probability of failure is substantially lower than that of all the other algorithms for any value of SNR. Our simulations referring to S3 have allowed us to assess how the performance of the considered algorithms is influenced by the overall number of tones. Some results referring to this scenario are shown in Fig. 9, which shows the dependence of P f on the number of tones at an average SNR = 10 dB. From this figure it is easily inferred that: 1) The CFH algorithm and the Alg-YA have a failure probability close to unity for L ≥ 5.  2) The CSFEC and the HCSFEC techniques are substantially more robust than all the other estimators, since their failure probability is lower than that of the other algorithms for any value of L.
3) The HCSFEC outperforms the CSFEC algorithm since the last estimator is biased (see Fig. 6). Finally, our analysis of S4 has allowed us to assess how different estimators perform in the presence of many closely spaced tones, having different strengths. The dependence of RMSE f , evaluated for the each of the normalized frequencies, on the tone spacing ∆F · N is shown in Fig. 10, whose results have been obtained under the assumption that SNR = 10 dB for the strongest tone (i.e., for the tone having the smallest frequency (see eq. (125) with k = 0)). The algorithms considered in this case are the CSFEC and the HCSFEC only, since the P f of the CFH algorithm and Alg-YA is close to unity (see Fig. 9). The results shown in Fig. 10 leads easily to the following conclusions: 1) The RMSE f characterizing the HCSFEC is lower than that achieved by the CSFEC for any value of tone separation and for any tone, since the CSFEC technique suffers from the error introduced by its frequency bias. 2) A larger RMSE f is observed for weaker tones at an arbitrary ∆F .
3) The RMSE f characterizing each tone reaches a floor for ∆F > 1.9/N .

B. Experimental Results about MIMO Radar Systems
A measurement campaign has been accomplished in the building of our institution to acquire a data set through a FMCW MIMO radar and a SFCW MIMO radar, both operating in the E-band. The first device is the TIDEP-01012 Cascade mmWave radar; it is manufactured by Texas Instrument Inc. [71] and is classified as a long range radar (LLR). Its main parameters are: a) chirp slope µ = 4 · This choice allows us to achieve the range resolution ∆R 2 = c/(2B) = 7.5 cm and the azimuthal resolution ∆θ 2 = 7.6.
All the measurements have been acquired through our radar devices in a large empty room (whose width, depth and height are 10 m, 8 m and 2.5 m, respectively) under the following experimental conditions: a) The two radar devices have been mounted on an horizontal wooden bar and have been lifted by a tripod at an height of roughly 1.60 m from ground (see Fig. 11); moreover, they have been activated alternately to avoid mutual interference.
b) The measurements acquired through the SFCW MIMO radar have been always pre-processed by the cancellation algorithm already available on this device; this algorithm exploits the measurements acquired from the first transmitted frame to cancel out unwanted received echoes. c) As already mentioned above, a small virtual ULA, consisting of N V = 16 virtual antennas, has been exploited in both devices. This has allowed us to reduce the far-field distance, which is proportional to the square of the overall size of the employed ULA.
In our experimental set-up, a pico-flexx camera manufactured by PMD Technologies Inc. [73] has been employed as a reference sensor; this device is based on a near-infrared vertical cavity surface emitting laser (VCSEL), and is able to provide a depth map or, equivalently, a three-dimensional 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 • ).
As far as the acquired measurements are concerned, it is important to point out that: a) In both the FMCW and SFCW radar systems, all the target ranges have been estimated with respect to the central virtual channel of their ULA.
b) The exact target positions have been acquired with respect to the centre of the pico-flexx camera. Therefore, in comparing these positions with their estimates computed on the basis of the radar measurements, the distance ∆ F P = 50 cm (∆ SP = 20 cm) between the FMCW (SFCW) radar and the camera was always kept into account. c) All our measurements have been processed in the MATLAB environment. A desktop computer equipped with a single i7 processor inside has been always used.
The numerical results illustrated in this paragraph refer to two specific static scenarios. The first scenario is characterized by a single detectable target, represented by a small metal disk 7 of size 5.5 cm. The target range R has been varied in the interval (1.0 m, 3.0 m), with a step of 0.5 m; in all cases, its azimuth θ has remained within the interval (−40 • , 40 • ). The exact positions considered for our target are listed in Table II for both the FMCW and the SFCW MIMO radar systems (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 5 (so that 1 ≤ L ≤ 5). As shown in Fig. 11, the targets are represented by small coins with a diameter of 2 cm, whose exact positions are listed in Table III  (the data referring to the i-th target are collected in the column identified by T i , with i = 1, 2, ..., 5); they been sequentially added in our scenario, starting from a single target and then increasing the number of targets up to 5. This has allowed us to assess how the performance of our estimation algorithms is influenced by the value of the parameter L. In processing all the acquired measurements, prior knowledge of L has been always assumed and the following values have been selected for the parameters of the CSFE (CSDE) and CSFEC (CSDEC) algorithms: Moreover, in analysing the data acquired in both scenarios, the following choices have been made: 1) The accuracy of range estimates has been assessed by evaluating the RMSĒ and the peak errorε where N m represents the overall number of available measurements.
2) 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 estimated by the CSFE (CSDE) over the ULA and associated with a target placed at approximately 8 θ = −14 • with respect to the center of the two radars is shown in Fig. 12. Since the distance d H between adjacent virtual channels is constant, the estimated phases exhibit a linear dependence on the index of the virtual channel (see eqs. (13) and (14)). 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 eq. (25)); this is confirmed by the results shown in Fig. 12 for both the FMCW and the SFCW radar systems. To assess the quality of the estimated phases, their RMSE 9ε ψ has been evaluated in both scenarios; in doing so, the linear fitting of the phases estimated over the whole ULA has been taken as a reference with respect to which each phase error has been computed.
The estimate of the target range generated by the CSFE and the CSDE in each of the N m = 10 distinct positions considered in the first scenario are listed in Table II; in the same table, the value ofε ψ computed for each position is also provided . These results have allowed us to compute the errorsε R (130) andε R (131), the mean ofε ψ (denotedε m,ψ and generated by taking the average of the N m values available forε ψ ) and the average computation time (CT); their values are summarized in Table IV for the same scenario.
These results and those listed in Table II have led us to the following conclusions: 1) Both the CSFE and CSDE are able to accurately estimate the range and the amplitudes of a single target.
2) The value ofε ψ characterizing the CSDE algorithm is lower when the target is closer to the radar, whereas the opposite occurs in the case of the CSFE. This result is mainly due to the fact: a) since the antenna array of the SFCW radar is more compact than that of the FMCW radar, the far-field condition of the former device is satisfied at shorter ranges than the second one: b) the SFCW (FMCW) radar is a SRR (LRR). In fact, for these reasons, the measurements provided by the SFCW radar allow to obtain more accurate results at a shorter distance. On the contrary, the higher transmit power radiated by the FMCW allows to achieve a smallerε ψ at larger ranges. 3) All the values ofε R andε R are comparable, reasonably low and in the order of the resolution of our devices. 4) The value ofε m,ψ achieved with the SFCW radar is slightly lower that obtained with the FMCW radar.
Methods   5) The CTs are in the order of few milliseconds. Let us focus on the second scenario now. The range estimates generated by the CSFEC, HCSFEC, CSDEC, HCSDEC, CFH and HASQE algorithms, and the Alg-YA in this scenario are listed in Table III. Note that the estimates provided for the CSFEC, HCSFEC, CSDEC and HCSDEC refer to 5 targets, whereas those provided for the CFH and HASQE algorithms, and for the Alg-YA to 4 targets only. This is due to the fact that the last three algorithms, unlike the first four, do not converge in the case of 5 targets. The errorsε R andε R , and the CT assessed in this case are listed in Table V. From these results it is easily inferred that: 1) The CSFEC (CSDEC) algorithm and the Alg-YA achieve better accuracy than the other algorithms. However, unlike the Alg-YA, the CSFEC (CSDEC) works properly in the presence of five targets. 2) In the considered scenario, the hybrid version of the CSFEC (CSDEC) algorithm does not provide a better accuracy than CSFEC (CSDEC); this means that estimation bias of the last algorithms does not play an important role in this case.
3) The estimated RMSEs and peak errors are in the order of the resolution of our devices, but are a little bit higher in the SFCW radar system. We believe that this is mainly due to the poorer estimation of the first target, since, in our specific experiment, the energy received from this target has been found to be lower than that coming from the other four targets. This problem is not found in the FMCW radar system, since this device, being a LRR, is able to achieve an higher SNR at its RX side for all the targets. Finally, we would like to stress that the robustness of CSFEC and CSDEC algorithms is related to the accuracy of the estimation and cancellation procedure they accomplish. This is exemplified by Fig. 13, where the initial amplitude spectrum of the signal received on the central virtual channel of the FMCW radar system in the second scenario and its residual, resulting from the cancellation of the spectral contributions due to the five detected targets, are shown.

V. CONCLUSIONS
In this manuscript, a novel algorithm for detecting and estimating a single tone has been derived; moreover, its has been shown how it can be exploited to estimate multiple tones through a serial cancellation procedure. The accuracy and robustness of the devised estimators has been assessed by means of extensive computer simulations and by analysing their application to range detection and phase estimation in commercial MIMO radars. Moreover, our results have evidenced that these estimators outperform all the other related estimators available in the technical literature in terms of probability of convergence and accuracy, in the presence of both a single tone and multiple tones. Future work concerns the application of the developed algorithms to other fields.

APPENDIX B A. Spectrum cancellation
In this paragraph, the expression of the vectorC 0 (Â l [r − 1],F l [r − 1]) appearing in the RHS of eq. (106) is derived. First of all, let us take into consideration the CSFE; in this case,C 0 (·, ·) is computed to cancel the contribution of the sequence (see eq. (91)) s n F l ,Ā l Ā l exp j2πnF l =Ā lw n l to the vector X 0 (42); here,w l exp(j2πF l ).
Since X 0 is the N 0 -th order DFT of the zero-padded vector x 0,ZP (43) (where the vector x 0 collects the elements of the complex sequence {x c,n ; n = 0, 1, ..., N − 1}), it is easy to prove that where Therefore, the identity holding for any q ∈ C, can be exploited for an efficient computation of all the elements of the vectorW (l) 0 . Similar considerations can be formulated for the SFE. However, in this case, eq. (140) is replaced by (see eq. (33)) s n F l ,C l C l exp j2πnF l +C * l exp −j2πnF l =C l (w l ) n +C * l (w * l ) n , wherew l is still defined by eq. (141). Consequently, eq. (142) is replaced bȳ whereW (l) 0 is the N 0 -dimensional vector defined above andW where Again, the identity (146) can be exploited for an efficient computation of all the elements of both the vectorsW Consequently, eq. (142) still holds, but has a different meaning, sinceW A similar procedure for leakage cancellation can be used also for the CSDEC algorithm. In this case, similarly to the CSFEC, the expression ofX lk,k is still defined by eq. (154); however, the quantityW

APPENDIX C
In this Appendix, the computational complexity, in terms of flops, is assessed for the SFE and the CSFE developed in Paragraph III-A. The overall computational cost of each these algorithms can be expressed as where C i,alg , C r,alg and N alg represent the cost of its initialization, that of each of its iterations and the overall number of iterations, respectively, and 'alg' denotes the algorithm these parameters refer to (SFE or CSFE). The general criteria adopted in estimating the computational cost of both C i,alg and C r,alg are the same as those illustrated in [74] and can be summarised as follows: • 4d − 2 flops are required to compute the inner product u T c v of a d × 1 complex column vector and a d × 1 real column vector; • 6d + 2(d − 1) flops are required to compute the inner product u T c v c of two d × 1 complex vectors; • d flops are required to find the maximum element of a vector v ∈ R 1×d ; • 4d 2 + 14d − 8 flops are required to compute an interpolation based on the elements of a complex vector v ∈ C 1×d ; The detailed expressions of C i,alg and C r,alg are provided below for both SFE and CSFE. SFE -The cost C i,SFE is evaluated as where: a) C X0 = 8N 0 log 2 N 0 is the contribution due to the computation of the vector X 0 (42); b) Cα = 4N 0 is the contribution due to the computation ofα on the basis of eq. (77); c) C Kp = 6 log 2 N + 151 is the contribution due to the evaluation of the quantities {K 1 (2α), K 2 (2α), K 3 (2α)} on the basis of eq. (62); d) C X k,α = 14N + 10 is the contribution due to computation ofX k,α on the basis of eq. (63); e) C∆ = 65 is the contribution due to the computation of the coefficients of the quadratic equation (69) and to the evaluation of its solution∆. The cost C r,SFE , instead, is evaluated as C r,SFE C g + CX + Cρ + CĈ + C X k,ρ + C Kp + C∆ where: a) C g = 15 is the contribution due to the computation of g(F ) on the basis of eq. (50); b) CX = 6N + 4 is the contribution due to the computation ofX(F ) on the basis of eq. (38); c) Cρ = 1 is due to the evaluation ofρ (i−1) on the basis of eq. (80); d) CĈ = 17 is the contribution due to the evaluation of the complex amplitudeĈ on the basis of eq. (49); e) C X k,ρ = 14N + 10 is the contribution due to the evaluation of the quantities {X k,ρ (i−1) } on the basis of eq. (63); f) the cost C Kp = 6 log 2 N + 151 is the contribution due to the evaluation of the quantities {K p (2ρ (i−1) )} on the basis of eq. (62) g) C∆ = 65 is the contribution due to the computation of the coefficients of the quadratic equation (69) and to the evaluation of its solution∆. It is worth mentioning that, as explained in Paragraph III-A, an interpolation technique can be used to compute the quantitiesX(F ) and {X k,ρ (i−1) }. In this case, CX = 4I 2 + 14I − 8 flops and C X k,ρ = 3(4I 2 + 14I − 8) flops are needed for the computation ofX(F ) and {X k,ρ (i−1) }, respectively, if a barycentric interpolation technique is used [69]. Substituting the terms appearing in the RHSs of eqs.
From the last expression, eq. (120) can be easily inferred. CSFE -The cost C i,CSFE is evaluated as where: a) C X0 and Cα have the same meaning as that already illustrated for the SFE; b) C X k,α = 33N + 15 is the contribution due to computation ofX k,α on the basis of eq. (63); b) C∆ = 29 is the contribution due to the computation of the coefficients of the quadratic equation (69) and to the evaluation of its solution∆. The cost C r,CSFE , instead, is evaluated as C r,CSFE CX + Cρ + CÂ + C X k,ρ + C∆ (166) where: a) CX and Cρ have the same meaning as that illustrated for the SFE; b) CÂ = 6N + 2 is the contribution due to the evaluation ofÂ on the basis of eq. (95); c) C X k,ρ = 33N + 15 is the contribution due to the evaluation of the quantities {X k,ρ (i−1) } on the basis of eq. (63); d) C∆ = 29 is the contribution due to the computation of the coefficients of the quadratic equation (69) and to the evaluation of its solution∆. Similarly as the SFE, an interpolation technique can be employed to compute the quantitiesX(F ) and {X k,ρ (i−1) }. In this case, CX = 4I 2 + 14I − 8 flops and C X k,ρ = 3(4I 2 + 14I − 8) flops are needed for the computation ofX(F ) andX k,ρ (i−1) , respectively, if a barycentric interpolation technique is used. Substituting the terms appearing in the RHSs of eqs. (165) and (166)