Deterministic Signal Processing Techniques for OFDM-Based Radar Sensing: An Overview

In this manuscript, we analyze the most relevant classes of deterministic signal processing methods currently available for the detection and the estimation of multiple targets in a joint communication and sensing system employing orthogonal frequency division multiplexing. Our objective is offering a fair comparison of the available technical options in terms of required computational complexity and accuracy in both range and Doppler estimation. Our numerical results, obtained in various scenarios, evidence that distinct algorithms can achieve a substantially different accuracy-complexity trade-off.

available methods can be divided in direct sensing methods and indirect estimation methods. The former methods extract target information from the received signal without compensating for the effect of the data payload conveyed by its useful component [5]; moreover, they typically exploit computationally intensive compressed sensing techniques (e.g., see [6], [7], [8], [9]). The latter methods, instead, rely on the knowledge of a preliminary estimate of the communication channel. The evaluation of this estimate requires compensating for the contribution of the transmitted channel symbols to the received signal (such symbols are always known at the receive side of a colocated radar; e.g., see [10]). In this manuscript, we focus on indirect methods only and investigate their use in a colocated 1 OFDM-based JCAS system equipped with a single TX and a single RX antenna (i.e., of single-input single-output (SISO) type). Moreover, we take into consideration different classes of indirect methods, namely: 1) discrete Fourier transform (DFT)-based or correlation-based methods; 2) subspace methods; 3) maximum likelihood (ML) based methods. It is worth pointing out that, even if various overviews on JCAS systems have appeared in the last three years [4], [5], [11], [12], [13], [14], [15], [16], [17], [18], [19], none of them provides a comparative analysis of the above mentioned methods for sensing in an OFDM-based radar system. Written with the aim of filling this gap, this manuscript offers a fair comparison in terms of accuracy and computational complexity of various algorithms belonging to the aforementioned classes and highlights their peculiarities and limitations.
The remaining part of this manuscript is organized as follows. In Section II, the processing accomplished in an OFDM-based radar system is summarized and the model of the signal feeding target detection and estimation algorithms is illustrated. Section III is devoted to the description of various relevant estimations methods, and to the assessment of their computational complexity. The analyzed techniques are compared in terms of accuracy and complexity in Section IV. Finally, some conclusions are offered in Section V.

II. SYSTEM AND SIGNAL MODELS
In this section, the processing accomplished in a SISO OFDM-based JCAS system is sketched; our main objective is providing the mathematical model of the transmitted signal and of the corresponding received signal in the presence of multiple targets. In our analysis, we focus on the transmission of a single frame, consisting of M consecutive OFDM symbols; the radio frequency (RF) waveform conveying the transmitted frame is (e.g., see [20, eq. (4)]) x RF (t) = ℜ exp (j2πf c t) where f c denotes the frequency of the local oscillator employed in the TX up-conversion and is the complex envelope of the transmitted signal conveying the mth OFDM symbol (with m = 0, 1, . . . , M − 1); here, q(t) is a windowing function (employed for pulse shaping), s m,n is the channel symbol carried by the nth subcarrier of the mth OFDM symbol (with n = 0, 1, . . ., N − 1), N is the overall number of subcarriers, f = 1/T is the subcarrier spacing, T is the OFDM symbol interval, T s ≜ T + T G is the overall duration of the OFDM symbol and T G is the cyclic prefix duration (also known as guard time [10]). In this manuscript, a rectangular windowing function is assumed, so that q(t) = 1 for t ∈ [−T G , T ] and q(t) = 0 elsewhere.
Let assume now that x RF (t) (1) is reflected by K point targets, and that the kth target (with k = 0, 1, .., K − 1) is located at the (initial) distance R k from the transmitter and moves at the radial velocity 2 v k with respect to it. It is not difficult to show that, in this case, the complex envelope of the signal received by the JCAS system is 3 where A k ≜ α k exp(−j2πf c τ k ) is the complex gain accounting for the attenuation, path loss and phase rotation due to the overall propagation delay τ k ≜ 2R k /c introduced by the kth target and α k is a positive parameter accounting for the attenuation and the path loss associated with the same target. Moreover, f D k = 2v k /λ is the Doppler shift due to the relative velocity of the kth point target with respect to the radar system, λ = c/f c is the wavelength of the radiated signal and w(t) is the complex additive Gaussian noise process affecting r(t). The signal r(t) (3) undergoes analog-to-digital conversion and DFT processing. An analytical model that describes the sequence generated by sampling r(t) in the mth OFDM symbol interval can be obtained as follows. First of all, the right-hand side (RHS) of (2) is substituted in that of (3). Then, extracting the portion associated with the mth OFDM symbol from the resulting expression and substituting t with t ′ = t − mT s yields s m,n · a n (−F r k ) ξ n f D k , t ′ ζ m,n f D k exp j2πn f t ′ + w(t ′ ), with q being either m or n if X is either D or r, respectively, and are the normalized target delay and the normalized Doppler frequency, 4 respectively, characterizing the kth target, and Sampling r m (t ′ ) (4) at the instant t ′ l ≜ lT /N (that is equivalent to sampling r(t) (3) at the instant t l ≜ t ′ l + mT s ), with l = 0, 1, . . . , N − 1, results in where D l (f D k ) ≜ exp(j2πf D k lT /N ) accounts for the so-called range migration effect due to the kth target Doppler (e.g. see [21]), is the Gaussian noise affecting r m,l (an additive white Gaussian noise (AWGN) model is assumed for the sequence {w l ; l = 0, 1, . . . , N − 1}). In the following, we assume that the target Dopplers {f D k } are sufficiently small and, more precisely, |f D k /f c | ≪ 1/(MN ) for any k, so that the factors ξ n,l (f D k ) and ζ m,n (f D k ) appearing in the RHS of (10) can be neglected; this leads to the simplified signal model The N signal samples acquired in the mth OFDM symbol interval undergo serial-to-parallel (S/P) conversion; this produces the N -dimensional vector r m ≜ [r m,0 , r m,1 , . . . , r m,N −1 ] T for which an order N DFT is computed. The nth element of the resulting DFT output vector can be expressed as where W m,n is the Gaussian noise affecting the nth subcarrier of the mth OFDM symbol and is a normalized frequency accounting for the target delay and the range migration due to its velocity. In the following we assume that N is large enough, so that the F ρ k ≈ F r k for any k, thus (13) can be simplified as Since the channel symbol s m,n is known at the receive side for any n and m, the estimatê is the noise sample affectingĤ m,n , in (16). It is worth pointing out that: 1) The parameters F r (6) and F D (7) satisfy the inequalities F r,min ⩽ F r ⩽ F r,max and F D,min ⩽ F D ⩽ F D,max , with F r,min = 0, F r,max = 1 and F D,min = −1/2, F D,max = 1/2, respectively.
2) an AWGN model is adopted for the noise samples {W m,n } (see (15)) and, consequently, for the noise samples {W m,n } (see (17)), since a phase shift keying (PSK) constellation is employed in our simulations (each element of the sequence {W m,n } is assumed to have zero mean and variance σ 2 W ).
3) Neglecting self-interference and range migration results in a signal model in which the target delay and Doppler frequency are decoupled parameters. This entails that, in principle, the values of these parameters can be evaluated separately through 1D frequency estimation techniques. However, only 2D frequency estimation techniques are considered in the following since, despite their higher computational effort than their 1D counterparts, they achieve better accuracy and do not require the use of a pairing method to associate the estimated delays and Doppler frequencies with each detected target. 4) Even if clutter plays an important role in radar sensing, its contribution to the received signal can be mitigated resorting to various techniques available in the technical literature (e.g., see [6]). For this reason, in our work, this contribution is always neglected.
From (16) it can be easily inferred that: a) the noisy samples {Ĥ m,n } of the two-dimensional (2D) channel response acquired over a single frame can be modelled as the superposition of multiple 2D complex exponentials with AWGN; b) target detection and estimation is tantamount to identifying the K complex exponentials forming the useful component of the 2D sequence {Ĥ m,n } and at estimating their parameters, respectively.

III. DETECTION AND ESTIMATION ALGORITHMS
In this subsection, various algorithms for the detection and estimation of multiple targets in an OFDM-based radar systems are illustrated and their computational complexity is analyzed by deriving the order of magnitude of the number of floating point operations (FLOPs) they require to process a single OFDM frame. The general criteria adopted in estimating the computational cost of the various algorithms are the same as those illustrated in [22] and [23,Appendix C]. These algorithms are divided in FFTbased techniques, subspace-based methods and ML-based techniques.
68874 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

A. FFT-BASED TECHNIQUES
Correlation-based and DFT-based methods have been developed in [24], [25], [26], [27], [28], [29], [30], [31], and [32]. In particular, matched filter (MF)-based techniques for the estimation of range and Doppler in a single or multi-target scenario have been investigated in [24] and [32]. Such techniques benefit from the prior knowledge of the received signal and are computationally efficient; however, the accuracy they achieve in radar imaging may be poor because of high sidelobes and leakage, especially in the presence of strong clutter around real targets. A serial cancellation technique for improving the overall accuracy of radar images has been developed in [33], whereas a reduced complexity method, based on the idea of splitting a 2D estimation problem (involving target range and Doppler) into a couple of onedimensional (1D) simpler sub-problems, has been illustrated in [29].
In the following we take into consideration the 2D periodogram method (dubbed 2D-FFT in the following) and two cancellation-based estimation algorithms, namely the complex single frequency delay estimation and cancellation (CSFDEC) algorithm [33] and the CLEAN algorithm [34], [35]. On the one hand, the first algorithm can be considered as a reference technique; on the other hand, the other two algorithms as methods able to efficiently mitigate the main problems of the first one, namely: 1) the limited accuracy due to the discretization of the grid selected in the search for the peaks of the periodogram; 2) the need to search for multiple local maxima, which can lead to missed target detection in the presence of closely spaced targets; 3) the significant impact that spectral leakage may have in the presence of multiple targets. In fact, both the CSFDEC and the CLEAN algorithms combine the serial cancellation of the spectral contribution of each detected target with leakage compensation and re-estimation techniques. Moreover, unlike the 2D-FFT algorithm, they do not need a prior knowledge of the overall number of targets. 5

1) TWO-DIMENSIONAL PERIODOGRAM METHOD
This method is based on the so-called range-Doppler map [26], i.e. on the function with l = 0, 1, . . . , M 0 − 1 and p = 0, 1, . . . , N 0 − 1; here, is the coefficient (l, p) of the order (M 0 , N 0 ) discrete symplectic Fourier transform (DSFT) of the 2D sequence {Ĥ m,n }. Moreover, and L D and L r are the oversampling factors adopted in the Doppler and range domain, respectively. The estimates of the normalized Doppler frequency F D and the normalized delay F r of a single target are evaluated asF and S X is the set of integers {0, 1, . . . , X − 1} for any positive integer X . Givenl andp, the target complex amplitude A is estimated asÂ In a multi-target scenario, multiple (sayK , whereK denotes a prior estimate of the number of targets K ) local maxima 6 should be expected in the range-Doppler map; in this case, the parameters of the target associated with each local maximum are evaluated according to (20), (21) and (25). The most computationally intensive task required by this method is represented by the evaluation of the above mentioned order (M 0 , N 0 ) DSFT (see (19)). Then, the cost for the search ofK local maxima in the range-Doppler map has to be added to the previous cost. For this reason, the overall computational cost of the 2D-FFT method is 2) CLEAN ALGORITHM The use of the CLEAN algorithm in radar systems has been proposed in [34], [35], and [36]. On the one hand, one of the earliest implementations of the CLEAN algorithm can be found in [34], where the cancellation capability of the CLEAN algorithm is adopted to reduce sidelobe-induced artifacts affecting the images generated by microwave systems that employ antenna arrays. On the other hand, a more recent CLEAN-based algorithm has been proposed in [35] and [36], where it is employed in the context of stepped frequency continuous wave (SFCW) radar technology. In that case, the algorithm also includes a technique for the compensation of the spectral leakage due to close targets. In this subsection, we show how this algorithm can be also employed for jointly estimating the range and velocity of multiple targets in the OFDM-based radar system described in the previous section. The CLEAN algorithm is based on the same cost function as the 2D-FFT method (see (18)), but, unlike it, makes use of an iterative target cancellation procedure. This means that, within each of its iterations, after detecting a new target and estimating its parameters, its contribution to the above mentioned cost function is cancelled; this results in a residual cost function, which is passed to the next iteration. More precisely, the processing executed by the CLEAN algorithm consists in an initialization step followed by an iterative procedure. In the initialization, we set the iteration index k to 0 and with m = 0, 1, . . . , M − 1 and n = 0, 1, . . . , N − 1. Then, in the kth iteration (with k = 0, 1, . . . ,K − 1, whereK denotes the overall number of detected targets), the four steps described below are carried out sequentially. 1) Computation of the cost function -The cost function is computed for l = 0, 1, . . . , M 0 − 1 and p = 0, 1, . . . , N 0 − 1; here, S k [l, p] is expressed by (19), where, however,Ĥ m,n is replaced byĤ m,n [k] (see (29) below).
2) Estimation of the parameters of a new target -A search for the global maximum over the set {J k [l, p]; l ∈ S M 0 , p ∈ S N 0 } (collecting M 0 N 0 values) is performed to detect a new target (i.e., the kth target); the value of the couple (l, p) corresponding to the global maximum is denoted (l k ,p k ). Then, the estimates of the normalized frequency F D k and of the normalized delay F r k of the kth target are evaluated aŝ , respectively (see (20) and (21)); whereas that of its complex amplitudeÂ k is evaluated according to (25) 3) Threshold test to identify false targets -If |Â k | < T CLEAN , where T CLEAN denotes a proper (positive) threshold, a false target is identified and the execution is stopped by moving to step 5); otherwise, we proceed with the next step.

4) Target cancellation -The new residual frequency responsê
is evaluated to cancel the contribution of the last detected target toĤ m,n [k] (with m = 0, 1, . . . , M − 1 and n = 0, 1, . . . , N − 1). Then, the iteration index k is increased by one and a new iteration is started (i.e., we go back to step 1)). 5) End -The final output provided by the CLEAN algorithm is represented by the set {(F D k ,F r k ,Â k ); k = 0, 1, . . . ,K − 1}, whereK represents the last value taken on by the iteration index k.
The serial cancellation procedure expressed by eq. (29) may suffer from error accumulation, since the effects of errors in the estimation of target parameters accumulate over successive iterations. This may result in: a) poor accuracy in the presence of multiple and/or closely spaced targets; b) the detection of false targets. These considerations motivate the use of a refinement procedure to be accomplished after the last iteration of the CLEAN algorithm. This procedure consists of of the parameters of theK targets detected by the CLEAN algorithm are evaluated as follows. First of all, we maximize, over a specific rectangular grid (consisting ofM 0Ñ0 distinct nodes), the function where with k = 0, 1, . . . ,K − 1. The above mentioned grid, that has a significant impact on the accuracy achieved by the CLEAN algorithm, has the following properties: 1) its center depends on bothF Table 1 it is shown how the nodes of the above mentioned grid are selected; in that Table, The nodes are grouped in the set where × denotes the Cartesian product between two sets, I It is important to point out that: 1) the evaluation ofĤ ); k = 0, 1, . . . ,K − 1} become available; 3) the value assigned to δ X (with X = D and r) allows to cover two adjacent bins of the evaluated DSFT.
It can be shown that the computational cost of the CLEAN algorithm with refinement is O(N CL ), where (see [36,Sec. III-E, eq. (43)]) 7 In the equations listed in Table 1 and in (36), the dependence ofF D,r on the target index k is not specified to ease notation.
is the contribution due a single iteration of the algorithm; note that the parameters (M 0 , N 0 ) and (M 0 ,Ñ 0 ) define the grid sizes for the initialization and for the refinement steps, respectively.

3) CSFDEC ALGORITHM
The CSFDEC algorithm, developed in [33], combines a single 2D tone estimator, named complex single frequency delay estimator (CSFDE), with a serial cancellation procedure, similar to that used by the CLEAN algorithm [35]. However, the CSFDEC algorithm performs target cancellation in the frequency domain, whereas the CLEAN algorithm accomplishes it on the 2D sequence {Ĥ m,n } extracted from the (time domain) received signal. The processing accomplished by the CSFDE algorithm consists in an initialization followed by an iterative procedure. In the initialization, the following quantities are computed: 1) The set of 13 (3,0) and (3,3)); here, is the order (M 0 , N 0 ) DSFT of the zero padded version and with m = 0, 1, . . . , M − 1 and n = 0, 1, . . . , N − 1 (the oversampling factors L D and L r are defined by (22) and (23), respectively).
2) The coarse estimatesF r,c = p (0) /N 0 of the normalized Doppler F D and normalized delay F r , respectively, that characterize the single target detected by the algorithm; here, 3) The initial estimatê of the complex amplitude A of the detected target; here, and whereȲ k 1 ,k 2 ≜Ȳ k 1 ,k 2 [l (0) ,p (0) ] for any k 1 and k 2 . 5) The initial estimateˆ (0) (ˆ (0) ) of ( ) aŝ with X = (with X = ). 6) The initial fine estimateŝ of F D and F r , respectively. This concludes the initialization. The iterative procedure is started by setting the iteration index i = 1. The ith iteration is fed by the estimatesF of F D , F r and A, respectively, and produces the new estimatesF r andÂ (i) of the same quantities, with i = 1, 2, . . . , N it (where N it is the overall number of iterations selected before starting the algorithm). The procedure adopted for the evaluation ofF r andÂ (i) consists of the two steps described below.
1) Estimation of the normalized Doppler and the normalized delay -The new estimatesˆ (i) andˆ (i) of˜ and˜ , respectively, are computed according to (49). In evaluating the coefficients b X and c X (with X = and ) on the basis of (45)-(48),Â =Â (i−1) is assumed; moreover, we set ) for any k 1 and k 2 , wherē An alternative to the use of the last formula for the evaluation ofȲ k 1 ,k 2 is represented by the use of a 2D interpolation method applied to a block of I D I r elements of the matrixȲ k 1 ,k 2 (where I D and I r denote the interpolation orders adopted for Doppler and range, respectively); the need of interpolation originates from the fact that, in general,F (20) and (21) (with a proper choice of the integers l and p), respectively. Givenˆ (i) andˆ (i) , the new estimateŝ are computed.
2) Estimation of the complex amplitude -The new esti-mateÂ (i) ofÂ is evaluated by means of (43); in doing so, the couple (F Then, before starting the next iteration, the index i is incremented by one. At the end of the last (i.e., of the N it th) iteration, the fine estimatesF D =F A =Â (N it ) of F D , F r and A, respectively, become available and the algorithm stops.
The CSFDE algorithm represents the core of the CSFDEC algorithm, which is used to recursively estimate the multiple tones forming the useful component of the 2D complex sequence {Ĥ m,n }, whose (m, n)th element is expressed by (16), with K ≥ 1 and, in general, unknown. The CSFDEC algorithm is initialized by: 1) running the CSFDE algorithm to compute the initial estimatesF of the parameters F D 0 , F r 0 and A 0 characterizing the first (i.e., the strongest) target; 2) setting the recursion index i to 1 andȲ (0) 0,0 =Ȳ 0,0 (see (39) with k 1 = k 2 = 0). Then, a recursive procedure is started. The ith recursion of this procedure is fed by the vectorsF where T CSFDEC is a proper threshold, the ith tone has to be considered as the last one; for this reason, the estimateK = i of K is generated and the algorithm stops after carrying out the next two steps. In the second step, N it iterations are executed to refine the estimate of the parameters of the new tone detected in the previous step. The processing accomplished in this step follows closely that described in the refinement part (i.e., in the second step) of the CSFDE. Therefore, in each iteration, a new estimate of the complex amplitude and of the two frequency residuals of the ith tone are computed. Finally, in the third step, each tone is re-estimated after cancelling the leakage due to all the other (i − 1) tones. This allows to progressively refine the amplitude, normalized Doppler frequency and normalized delay of each tone, thus generating the final estimates of all the detected tones. Note that, in principle, this re-estimation procedure can be repeated multiple (say, N REF ) times.
Finally, it is worth noting that the refinement procedure of the CSFDEC algorithm does not need, unlike that of the CLEAN algorithm [35], the definition of a search grid for the optimization of a cost function. In fact, the frequency residuals are computed by the CSFDEC on the basis of a closed from expression (see (49)), whose use requires only the knowledge of the spectral coefficients collected in the matrices {Ȳ k 1 ,k 2 }. It can be shown that the computational cost of the CSFDEC algorithm, when the frequency residuals are evaluated on the basis of (49), is O(N CSFDEC ), where (see [33, Sec. III-C, eq. (58)]) (55)
In particular, algorithms based on the Multiple Signal Classification (MUSIC) technique have been analyzed in [9], [20], [37], and [39], whereas the use of the 2D estimation of signal parameters via rotational invariant technique (ESPRIT) for the estimation of the delay and Doppler of multiple targets has been studied in [38]. Various results illustrated in the above mentioned manuscripts lead to the conclusion that MUSIC-based algorithms can outperform 2D FFT-based methods for joint range-velocity estimation at the price, however, of a significantly larger computational complexity [37]. A lower complexity version of the MUSIC technique, called auto-paired method, has been proposed in [39], whereas an iterative method for improving the robustness of the MUSIC algorithm has been developed in [9]. In the following we concentrate on the 2D-MUSIC algorithm only, since it performs similarly as the 2D-ESPRIT at a comparable computational cost [40], [41]. This algorithm is based on the search ofK local maxima in the 2D-MUSIC spectrum (also known as pseudo spectrum), whose computation requires the identification of the so called noise subspace. In practice, the algorithm needs a prior estimate (denotedK ) of K and consists of the following steps: 1  (16)).
2) The pseudo-spectrum is  (20) and (21), respectively) and Q n is (MN ) × (MN −K ) matrix, whose columns represent the (MN −K ) eigenvectors of R (56) associated with noise. 8 3) The estimates of the normalized delay and of the normalized Doppler frequency are evaluated by finding the global maximum of the pseudo-spectrum (57) in the case of a single target or itsK highest local peaks in the case of multiple targets. Note that, in the last case, all the local maxima of the cost function can be really identified if the spacing between all the targets is greater than the grid step size; this problem can be mitigated by using a finer grid (i.e., a smaller step size) at the price, however, of a larger computational complexity.
The overall 2D-MUSIC complexity is O(N MUSIC ), where Here, we have that: a) C R = 10(MN ) 2 + 2MN − 2 is the cost due to the computation of the covariance matrix R (56); b) C e ≈ (MN ) 3 is the contribution due to computation of the eigenvalues of the matrix R (56); c) C P = M 0 N 0 (8(M N ) 2 + M N ) is the contribution due to the evaluation of the pseudospectrum (see (57)).

C. MAXIMUM LIKELIHOOD-BASED TECHNIQUES
Maximum likelihood-based algorithms are able to achieve an accuracy comparable to that of subspace methods at the price of an increase in computational complexity [41], [42], [43], [44]. Recent research contributions to this field concern: 1) the exploitation of alternating maximization approach to mitigate the computational complexity of ML estimation [41]; 2) the derivation of an iterative non-linear kernel least mean square (KLMS) based estimation technique for the estimation of target range [42]; 3) the development of an ML method, based on a kinematic model of detected targets, for estimating target speed [44].
In the following, we first take into consideration the approximate ML methods recently proposed in [41] and show how they can be adapted to our signal model 9 (16); the resulting algorithms, are dubbed modified Zhang ML (MZML) and modified alternating projection ML (MAP-ML). Both these algorithms are iterative; however, the first one maximizes a 2D cost function, whereas the second one exploits the method of alternating projections to turn a 2D optimization problem into a couple of simpler 1D sub-problems in order to mitigate the overall computational effort. Finally, we take into consideration the modified Wax & Leshem (MWL) algorithm developed in [35] and [36] for joint range and azimuth estimation of multiple targets in FMCW radar systems, and the expectation maximization (EM) algorithm developed in [36] and [41]. In particular, the last algorithm is exploited as follows. It is initialized be means of the 2D-FFT technique. Then, in each of its iterations, it accomplishes leakage compensation in the same way as the MZML algorithm. The resulting algorithm is denoted Modified Zhang EM (MZEM) in the following.

1) MODIFIED ZHANG MAXIMUM LIKELIHOOD ALGORITHM
This algorithm operates in an iterative fashion; in each of its iterations, the estimates of the detected targets are refined. It is initialized by: 1) setting the iteration index i to 1; 2) selecting the initial estimate of the normalized Doppler frequency F D k and that of the normalized delay F r k asF (0) k ], respectively (see (20) and (21) that collect the estimates of the normalized Doppler frequency and the normalized delay, respectively, of theK tones estimated in the previous (i − 1) iterations, and produces the new estimatesF r . Moreover, in the ith iteration, two distinct steps, that are repeated sequentially for each target (i.e., for k = 0, 1, . . . ,K − 1), are executed; the description of these steps is provided below for the kth target.
1) Computation of the cost function -In this step, the cost function 9 Note that, in our model, unlike the one considered in [41], inter-pulse and inter-subcarrier Doppler effects are neglected. VOLUME 11, 2023 is evaluated for (F D k ,F r k ) ∈ I is the (MN ) × (MN ) orthogonal projection matrix, is a (MN ) ×K matrix. 11 The uth element of the last matrix, namely M   [35] in Section III-A2 (see (33)-(35), in Table 1).
2) Target parameter estimation -In this step, the esti-matesF (i) D k andF (i) r k of F D k and F r k are computed as Then, the new estimateÂ (i) k of the complex amplitude A k is evaluated aŝ k (see below) on the trial variablesF D k andF r k is not always explicitly shown for simplicity. For the same reason, the dependence ofM The overall computational cost of MZML algorithm can be expressed as where C 0,MZML is the cost of its initialization (i.e., the same cost as the 2D-FFT method; see (26)), whereas represents the cost of a single iteration; here, we have that: 1) is the cost due to the evaluation of the projection matrix P

2) MODIFIED ALTERNATING PROJECTION MAXIMUM LIKELIHOOD ALGORITHM
The MAP-ML algorithm is initialized exactly in the same way as the MZML algorithm (see [41]) and employs the cost function J in the estimation of the parameters of the kth target; however, the method of alternating projection is exploited for the maximization of that cost function in order to replace the 2D optimization problem (64) with two 1D optimization problems. In practice, in its ith iteration, the frequency estimatesF (i) D k andF (i) r k appearing in the left-hand side of (64) are evaluated aŝ respectively (the sets I  (69) and (70) requires computing the matrices M (63)), respectively; 2) in the first iteration (i.e., for i = 1),F (0) D k is computed according to (20), with l =l The overall computational cost C MAP−ML of the MAP-ML algorithm can be expressed as where C 0,MAP−ML is the contribution due to its initialization (equal to C 0,MZML ; see (66)), whereas is the contribution due to each of its iterations. Moreover, in the last formula, we have that: 1) C P r ≈ 8N 0K 2 M 2 N 2 and C P D ≈ 8M 0K 2 M 2 N 2 are the costs due to the evaluation of the projection matrix P (i) k (·, ·) (62) in the first and second 1D optimization, respectively; 2) C J r ≈ 8N 0K M 2 N 2 and C J D ≈ 8M 0K M 2 N 2 are the costs due to the evaluation of the function J 3

) MWL ALGORITHM
Similarly as the CLEAN algorithm [35], the MWL algorithm operates in an iterative fashion and, in each of its iterations, estimates the parameters of a new target, and performs cancellation and leakage compensation in the time domain. However, unlike the CLEAN algorithm, the MWL algorithm requires solving 1D optimization problems only. Moreover, it can achieve similar and even better accuracy than the CLEAN algorithm with a smaller computational effort [36]. The MWL algorithm is initialized by setting the iteration index k to 0 andĤ m,n [0] =Ĥ m,n (see (16)), with m = 0, 1, . . . , M − 1 and n = 0, 1, . . . , N − 1. Then, its iterations are started; the processing accomplished in the kth iteration evolves through the four steps described below.
1) Coarse estimation of the Doppler of a new target -In this step, the coarse estimateF D k of the Doppler frequency F D k of a new (namely, of the kth) target is computed by solving the 1D optimization problem where F D [l] is defined by (20), a(F D [l]) is an M -dimensional column vector whose mth element (with m = 0, 1, . . whereĤ m,n [k] is evaluated on the basis of (29) if k > 0, with m = 0, 1, . . . , M −1 and m ′ = 0, 1, . . . , M −1. Givenl k (74), the coarse estimate of F D k is computed asF D k = F D [l k ] (see (20)).
2) Estimation of target delay -In this step, an estimateF r k of the normalized delay F r k characterizing the kth target is evaluated by solving another 1D optimization problem. This requires: a) Computing the N -dimensional column vector wherê where F r [p] is defined by (21), and a(−F r [p]) is an N -dimensional column vector whose n-element (with n = 0, 1, . . . , N − 1) is a n (−F r [p]) (see (5)). Givenp k (78), the final estimateF r k of the target delay is evaluated according to (21) with p =p k .
3) Fine estimation of target Doppler -In this step, the fine estimateF D k of the normalized Doppler F D k characterizing the kth target is evaluated by solving the last 1D optimization problem and, in particular, asF D k = F D [l k ] (see (20)), wherê B(F r k ) is an M -dimensional row vector, whose mth element is defined aŝ andĈ 4) Estimation of target complex amplitude -In this step, the estimateÂ of the complex amplitude A k is evaluated. Similarly as the CLEAN algorithm [35], at the end of the last step a false target is detected if |Â k | < T MWL , where T MWL denotes a proper (positive) threshold. When this occurs, the execution is stopped; otherwise, a new iteration is started going back to step 1).
An iterative procedure can be employed to refine the target estimates generated by the MWL algorithm. Similarly as the CLEAN algorithm, this procedure is based on: 1) estimating again the parameters of each target after removing the spectral contribution of the other (K − 1) targets, i.e. their leakage; 2) shrinking the search grid as iterations evolve (see (33)-(35), in Table 1). In the ith iteration of the refinement (with i = 1, 2, . . . , N REF ), the finer estimates {(F  VOLUME 11, 2023 of the normalized Doppler frequency for the kth target; here, I (i) D (M 0 ) denotes the set ofM 0 trial values selected forF D k in the ith iteration and R (i) [k] is an M ×M matrix whose element (m, m ′ ) is still defined by (75), where, however,Ĥ m,n [k] is replaced byȞ Then, we computê where I (i) r (Ñ 0 ) denotes the set ofÑ 0 trial values selected for F r k in the ith iteration, a(−F r k ) is an N -dimensional column vector whose nth element is a n (−F r k ) (see (5)),v (i) k plays the same role asv k in (78), but its nth element is defined as (see (77) where H represents the cost due to a single iteration of the algorithm; note that the parameters (M 0 , N 0 ) and (M 0 ,Ñ 0 ) define the grid size for the initialization and for the refinement step, respectively.

4) EM-BASED ALGORITHM
The EM algorithm [45] can be employed jointly with each of the algorithms described above to refine its estimates of target parameters [36], [41]. For this reason, we can assume that, in general, the EM algorithm is fed by theK - , collecting the initial estimates of the normalized Doppler, normalized delay and the complex amplitude of theK detected targets.
The EM algorithm operates in an iterative fashion; in each of its iterations, it executes an expectation step (E-step) followed by a maximization step (M-step). In our description of such steps, we focus on the ith iteration (with i = 1, . . . , N REF , where N REF denotes the overall number of iterations) and consider the kth target (with k = 0, 1, . . . ,K − 1). At the beginning of this iteration, the estimates (F ) are available for the normalized Doppler, the normalized delay and the complex amplitude, respectively, of the considered target ). The two steps accomplished within the considered iteration are described below.
1) E step -In this step, the cost function 2) M step -The new (and, hopefully, finer) estimatesF r k of F D k and F r k , respectively, are computed as In the equation above, the term J (i) EM k (·, ·) evaluated through (90) represents the cost function computed over a specific rectangular grid, defined by the trial values (F D k ,F r k ), whose center depends on bothF (i−1) D k andF (i−1) r k , and whose step sizes gets smaller as i increases. The grid employed in this case is generated according to the same criteria illustrated for the refinement procedure developed for the CLEAN algorithm (see (33)- (35), in Table 1 and the comments related to them).  The overall computational cost of a single iteration of the EM algorithm in the presence ofK targets can be expressed in a similar way as [36, Appendix C], i.e. as EM k (·, ·) (90), whereas C opt = 4M 0 N 0 is the cost of its optimization; 4) C A = 8MN is the cost due to the computation ofÂ (i) k on the basis of (93). In our work, as already mentioned in Subsection III-C, the EM algorithm has been employed to refine the estimates generated by the 2D-FFT algorithm. The overall computational complexity of the resulting algorithm, called MZEM, is O(N MZEM ), where N MZEM = N 2D−FFT + N EM (see (26) and (94)).
The computational complexity orders of the estimation algorithms described above and considered in our simulations are listed in Table 2.

IV. NUMERICAL RESULTS
The estimation algorithms described in the previous section have been compared, in terms of accuracy and computational effort, in four distinct scenarios for different values of the signal-to-noise ratio where σ 2 W represents the variance of each element of the complex noise sequence {W m (n)} (see (17)). The first scenario (S1) is characterized by a single target (i.e., by K = 1), whose and respectively, and its amplitude A k is set to 1 for any k; here,   square error whereX k [t] denotes the estimate of the parameter X k evaluated for the kth target in the tth Monte Carlo run, N mc is the overall number of runs and X = r (X = v) if target range (target velocity) is considered. Moreover, the following choices have been made: 1) in the evaluation of the RMSE r , the vector collecting target ranges and the one collecting their estimates have been organized according to an ascending order (from minimum to maximum range); 2) in the evaluation of the RMSE v , the vector collecting target velocities and the one collecting their estimates are sorted in the same way as the vectors referring to target ranges; 3) an SNR belonging to the interval [−20, 20] dB has been considered; 4) the number of targets (K ) has been always assumed to be known (so that events of missed detection are avoided). Note that the knowledge of K does not prevent all the algorithms from identifying false targets; in our work, unwanted detections contribute to the evaluation of the RMSE.
In the third scenario (S3), the range and velocity of the targets are computed according to the same strategy adopted in S2, but the SNR is fixed to 0 dB and the value of K ranges from 1 to 10. In this case, we focus on the computational effort required by all the algorithms considered in S1 and S2 and, in particular, we assess both their computation time (CT) and the overall number of FLOPs they require. 12 The fourth scenario (S4) is characterized by the same number of targets as S2 (i.e., K = 4) and by an SNR equal to 0 dB. Moreover, target ranges and speeds are computed according to (96) and (97), but smaller bin spacings (more precisely,R = 1 andv = 1) are assumed. For this reason, in this scenario, spectral leakage may substantially affect target estimation; this can be easily inferred from Fig. 1, where  the range-Doppler map (or ambiguity function) is shown for the considered case. In our analysis of S4, we assess the convergence speed of six iterative algorithms (namely, the CSFDEC, CLEAN, MWL, MZML, MAP-ML and MZEM algorithms) and, in particular, we analyze how their accuracy changes as the overall number of their iterations 13 ranges from 1 to 5.
Our interest in the four scenarios defined above can motivated as follows. The numerical results obtained in the first two scenarios show how the considered algorithms perform in the presence of a single tone and of multiple (but adequately spaced 14 ) tones, respectively, whereas those obtained in the third scenario allow us to assess the trend of their computational requirements when the overall number of targets increases. Finally, the fourth scenario sheds some light on the trade-off between estimation accuracy and computational effort of the considered iterative algorithms in the presence of closely spaced targets.
In our simulations, the following choices have been made. First of all, prior knowledge about K has been assumed and, unless differently stated, the following values have been selected for the parameters of the OFDM modulation 15 : 1) overall number of subcarriers N = 32; 2) overall number of OFDM symbols/frame M = 32; 3) subcarrier spacing f = 250 kHz; 4) cyclic prefix duration T G = 0.25 T = 1µs (consequently, the OFDM symbol duration is T s = 1/ f + T G = 5 µs); 5) carrier frequency f c = 79 GHz. These values entail that R bin = 18.7 m and v bin = 12 m/s. Secondly, the following values have been selected for the parameters of the considered algorithms: 1) the oversampling factor L D = 16 and L r = 16 have been chosen for Doppler and range estimation, respectively, in both the 2D-FFT and 13 In the case of the CSFDEC algorithm, each iteration corresponds to the execution of the re-estimation procedure described in Subsection III-A3. 14 This means that the spectral leakage affecting each tone and originating from all the other tones is limited; therefore, the frequency estimation error due to this phenomenon is negligible. 15 The choices we made for the following parameters have been dictated by the technical literature on OFDM-based JCAS systems (e.g., see [41]).   CSFDEC algorithms (so that M 0 = ML D = 512 and N 0 = NL r = 512; see (22) and (23), respectively); 2) the 2D-FFT method has been used to compute the initial estimates VOLUME 11, 2023 considered SNR range; 3) the accuracy of each algorithm attains the CRLB above its SNR threshold, but, at high SNRs, may reach a floor. The last phenomenon, observed for the 2D-FFT and 2D-MUSIC algorithms, is due to their limited accuracy 19 (which depends on the overall number of FFT bins and on the overall number of trial points employed in the computation of the steering vectors, respectively). The other algorithms, instead, take advantage of their refinement cycles, which improve the accuracy of their final estimates.
In addition, in analyzing the results shown in Fig. 2, readers should keep in mind that: 1) the computational complexity of the CSFDEC, CLEAN and MZEM algorithms is approximately 13, 2.3 and 1.35 times higher, respectively, than that of the 2D-FFT, whereas that of the MWL algorithm is very close to it; 2) the complexity of the 2D-MUSIC, MAP-ML and MZML algorithms is 352, 183 and 819 times higher, respectively, than that of the 2D-FFT.
Most of the considerations illustrated above for Fig. 2 also apply to the results shown in Fig. 3 for N = 32 and in Fig. 4 for N = 256. In both cases S2 is considered, but the results illustrated in Fig. 4 concern range estimation only and refer to a subset of the considered algorithms (2D-MUSIC, MZML and MAP-ML are ignored because of their huge memory requirements). Note that: 1) Unlike S1, an SNR threshold is visible in Fig. 3 for all the algorithms, mainly because of the presence of spectral leakage.
2) Independently of the considered algorithm, an higher number of subcarriers allows to achieve a lower RMSE r and reduce the SNR threshold.
3) In analyzing the results shown in Fig. 3, it should be kept into account that the computational complexity of the CSFDEC, CLEAN, MWL and MZEM algorithms is approximately 34, 177, 2.7 and 4.7 times higher, respectively, than that of the 2D-FFT. Moreover, the complexity of the 19 In the case of the 2D-FFT algorithm, this phenomenon can be mitigated by interpolating its cost function to improve the accuracy in peak detection. 2D-MUSIC, MAP-ML and MZML algorithms is 7532, 1518 and 7070 times higher, respectively, than that of the 2D-FFT.
The CT 20 and the overall number of FLOPs assessed in S3 are shown in Fig. 5 and Fig. 6, respectively. From these results and from those illustrated in S1 and S2, it can be easily inferred that: 1) the gap, in terms of CT, between the CLEAN, MAP-ML and 2D-MUSIC algorithms is small; 2) the MAP-ML algorithm takes advantage of the alternating projections method to reduce the computational effort of the MZML algorithm; 3) the CT of the MWL algorithm is the closest to that of the 2D-FFT algorithm, but, as shown in Fig. 3, the former algorithm achieves a poorer estimation accuracy than the latter one for SNR ∈ [−15, 3] dB; 4) the CSFDEC algorithm stays in the middle between high and low complexity algorithms, but is able to achieve excellent estimation accuracy, as evidenced by the RMSE results illustrated for the previous scenarios.
Some numerical results 21 referring to S4 are shown in Fig. 7 for RMSE r and RMSE v and in Fig. 8 for the CT. These results deserve the following comments: 1) the estimation accuracy of most of the considered algorithms does not improve after the first three iterations; 2) the CLEAN, MZEM and CSFDEC algorithms exhibit similar CTs if three iterations are carried out for each of them; 3) the MWL algorithm requires a lower CT, but generates poorer range and velocity estimates.

V. CONCLUSION
In this manuscript, eight different algorithms for detecting multiple targets, and for jointly estimating their range and velocity in a SISO OFDM-based JCAS system have been described. All belong to the class of indirect methods and are deterministic; moreover, three of them (namely, the CLEAN, MWL and CSFDEC algorithms) exploit an iterative cancellation procedure in the estimation of target parameters. Our numerical results evidence that, in the presence of a single target or of multiple well spaced targets, all the considered algorithms achieve reasonable accuracy. However, their behavior substantially changes in the presence of multiple closely spaced targets; in such conditions, our attention has focused on how accuracy and computational effort are influenced by the overall number of a) targets and b) refinement cycles employed to improve the quality of the final estimates of target parameters. The following conclusions can be formulated in the light of our simulation results: 1) the CLEAN and CSFDEC algorithms require a limited complexity and achieve excellent performance thanks to the use of a target cancellation procedure; 2) the MZML and MAP-ML algorithms achieve better accuracy than subspace-based method (i.e., than the 2D MUSIC algorithm) at the price of similar (and really high) computational complexities; 3) the MWL algorithm achieves a better accuracy-complexity trade-off than the MAP-ML algorithm. Therefore, in future OFDMbased JCAS system, the selection of a target detection and estimation algorithm requires a careful assessment of the pros and cons characterizing the various options available in the technical literature. Our ongoing work concerns the possible extensions of the considered algorithms to MIMO OFDMbased radars.