Joint Range and Doppler Estimation Using Spectrally Efficient FDM
Full text
Joint Range and Doppler Estimation Using Spectrally Efficient FDM Michele Mirabella‡,Member, IEEE, Pasquale Di Viesti‡,Member, IEEE, Christos Masouros†,Fellow, IEEE and Giorgio M. Vitetta‡,Senior Member, IEEE Abstract This paper explores the use of data modulated by a spectrally efficient frequency division multiplexing (SEFDM) waveform for sensing. We first show that, if the presence of a cyclic prefix is assumed in the transmitted signal, the problem of multiple target sensing is tantamount to the detection and estimation of an unknown number of complex two-dimensional complex tones. Then, a novel iterative estimation method, based on a maximum likelihood approach, is developed to solve the last problem. Our simulation results evidence that SEFDM represents a valid technical option over static or slowly varying channels, and that the proposed method achieves a better accuracy-complexity trade-off than other estimation techniques available in the technical literature. In particular, our numerical results show that, in various scenarios, the proposed method achieves a 15% –50% improvement in estimation accuracy with respect to other techniques with a limited computational complexity. Index terms— Cyclic prefix, Maximum Likelihood Estimation, Orthogonal Frequency Division Multiplexing, Spectral Efficiency, Radar Processing I. INTRODUCTION Wireless communication and radar sensing have been advancing independently for many years, despite sharing various similarities in terms of both signal processing and system architecture. In the last few years, substantial research efforts have been devoted to the design of wireless systems able to perform communication and radar functions jointly. The interest in such a class of systems, that accomplish integrated sensing and communication (ISAC), has been motivated by the advantages they offer in terms of device size, power consumption, cost and spectral efficiency (SE) with respect to traditional wireless systems in various applications [1]. A crucial challenge in ISAC is the design of waveforms allowing to achieve good performance in both communication and sensing under specific spectral and complexity constraints [2]. In this context, one of the available (and underexplored) options is represented by spectrally efficient frequency division multiplexing (SEFDM) [3], a modulation format closely related to orthogonal frequency division multiplexing (OFDM). In fact, similarly to OFDM, the SEFDM format puts into practice the frequency multiplexing concept by employing multiple subcarriers; however, unlike OFDM, it improves SE by giving up the orthogonality constraint, i.e., by increasing the spectral overlap among distinct subcarriers, to allow more subcarriers in a given bandwidth. Unfortunately, this advantage is obtained at the price of the unavoidable introduction of inter-carrier interference (ICI) at the receive (RX) side and, consequently, of a higher data detection complexity. This has motivated the interest in developing computationally efficient detection and decoding techniques for mitigating the impact of ICI in SEFDM communication systems [4]–[8]. Note that all these contributions address the problem of data detection in SEFDM from a communication perspective; in this manuscript, instead, we focus on the use of SEFDM signaling for wireless sensing, i.e., for estimating the range and radial speed of multiple targets in place of the data symbols conveyed by the radiated waveform. While the use of the orthogonal counterpart of SEFDM (namely, OFDM) has received significant attention in ISAC (e.g., see [1], [2], [9]–[16]), no attention has been paid to the adoption of SEFDM for similar purposes. Among these contributions, the work in [16] offers an overview of single-input single-output (SISO) OFDM radar techniques for target estimation, providing a unified perspective on existing methods. In contrast, the present work explores SEFDM from a sensing perspective, with the aim of partially filling this knowledge gap. In particular, we analyze its potential as a radar waveform and assess the inherent trade-off between sensing performance and communication efficiency according to the amount of compression introduced in subcarrier spacing. It is also worth mentioning that target detection and estimation algorithms developed for passive sensing in OFDMbased systems can be divided into direct and indirect methods. In general, direct sensing methods estimate radar target parameters directly from the RX signal without compensating for the communication payload. These approaches often rely on computationally intensive compressed sensing techniques [1]. Conversely, indirect sensing methods require an initial estimation of the communication channel, which is then used to subtract the contribution of the transmitted data symbols (e.g., see [9, Eq. (20)]). Indirect methods can be further categorized into: 1) discrete Fourier transform (DFT)-based or correlation-based ‡University of Modena and Reggio Emilia, Dept. of Engineering “Enzo Ferrari”, Via P. Vivarelli 10/1, 41125 Modena (Italy) and with Consorzio Nazionale Interuniversitario per le Telecomunicazioni (CNIT), †University College London, U.K., email: [email protected], pasquale.di[email protected], [email protected], [email protected]. This work has been supported by the European Union under the Italian National Recovery and Resilience Plan (PNRR). Specifically, it has been carried out within the “Telecommunications of the Future” partnership (PE00000001 - program “RESTART”), CUP E63C22002040007 - D.D. n.1549 of 11/10/2022, and within the initiatives of Mission 4 Component 2, Investment 1.4 (D.D. 1033 17/06/2022, CN00000023). Furthermore, the project has involved the MOST – Sustainable Mobility National Research Center.
2 methods [10]; 2) subspace methods [11], [12]; 3) maximum likelihood (ML)-based methods [13], [14]. A common feature across the aforementioned techniques is the use of an initial coarse estimate (typically obtained through the periodogram method applied to range-Doppler maps [15]), whereas the main differences among them concern their refinement stage, where different signal processing techniques are used to improve the target initial estimates. This paper is motivated by our interest in investigating the use of SEFDM for radio sensing. For this reason, we explore its application to the detection and estimation of multiple targets, both in range and Doppler, within an ISAC framework. The scope of this manuscript is threefold. First, we develop a simplified model for the received signal in a co-located singleinput single-output (SISO) SEFDM transceiver. In doing so, we account for the presence of a cyclic prefix (CP), typically overlooked in SEFDM studies, in the transmitted signal. Note that incorporating a CP eliminates inter-symbol interference (ISI), significantly easing channel estimation at the price of reduced SE. Second, we exploit the above-mentioned signal model and develop a novel approximate ML-based method for the detection and estimation of multiple targets. This method, dubbed Newton-based multiple cisoid refiner (NMCR), combines a coarse initialization based on the periodogram method with a Newton-based iterative refinement strategy that jointly optimizes target parameters (for this reason, it belongs to the class of indirect sensing techniques). Thirdly, we assess the accuracy and computational requirements of the NMCR method in multiple heterogeneous scenarios and we compare it with other methods available in the technical literature. Our numerical results lead to the conclusion that, in various scenarios, the NMCR can achieve 15% –50% improvement in estimation accuracy with respect to other techniques with reasonable complexity. The remaining part of this manuscript is organized as follows. In Section II, the processing accomplished in a SEFDM-based radar system is described, a simplified model is developed for the received signal feeding the NMCR algorithm and, based on this model, a simple method for the estimation of the channel matrix is proposed. Section III is devoted to the derivation of the NMCR algorithm and to the assessment of its computational complexity. The NMCR algorithm is compared, in terms of accuracy and complexity, with other estimation algorithms in Section IV. Finally, some conclusions are offered in Section V. Notation: Throughout this manuscript, the following notation is adopted: 1) (·)∗and (·)Hdenote the complex conjugate and the complex conjugate transpose (Hermitian operator), respectively; 2) modB[·]indicates the modulo Boperator (where Bis a positive integer); 3) ∗denotes the linear convolution operator between two signals or functions; 4) ℜ{x}and ℑ{x}indicate the real part and imaginary part, respectively, of the complex variable x; 5) The symbols ⊙and ⊘represent the Hadamard product and Hadamard division operators, respectively; 6) The symbol ⊛denotes the Khatri-Rao product operator; 7) ΞVis the unitary DFT matrix of order V, whose element (p, q)is exp(−j2πpq/V )/√V; 8) X≜[xm,n]defines a matrix Xof proper size and xm,n denotes the element appearing on its mth row and nth column; 9) y= vec(Y)defines an (MN)-dimensional column vector resulting from the ordered concatenation of the columns of the M×Nmatrix Y; 10) INis the identity matrix of order N; 11) diag(xN)generates a diagonal N×Nmatrix having the elements of the N-dimensional vector xNon its main diagonal; 12) Y†≜(YHY)−1YHis the Moore-Penrose pseudo-inverse of matrix Y; 13) Tr{Y}defines the trace of the matrix Y. II. SYSTEM AND SIGNAL MODELS In this section, the overall processing required for the sensing task in a SISO SEFDM-based ISAC system is briefly described; in doing so, emphasis is put on the processing at the RX side (which coincides with the transmit, TX, one in a co-located system). For better clarity, our description is divided into two parts. In the first part, the mathematical model of the received signal in the presence of multiple targets is developed and some essential assumptions on which it relies are illustrated. In the second part, instead, we show how channel estimation can be performed at the TX side of the system; this is then used to generate the signal to be employed for target detection and estimation. A. Received signal in the presence of multiple targets In the following, we take into consideration the transmission of a single SEFDM frame, consisting of Mconsecutive SEFDM symbols, over a slowly varying wireless channel. We assume that the SEFDM frame incorporates both pilot tones (for channel estimation and synchronization in digital communications) and information data to be sent to a single or multiple receivers at different locations. However, since the receiver is assumed to be co-located with the transmitter (i.e., monostatic sensing is employed), full knowledge of the structure and content of the whole frame, and of the transmission frequency is available at the RX side; this information is exploited by the SEFDM-based ISAC receiver for sensing purposes only, as shown in Fig. 1. Note that, in the following, we focus on a SISO radar architecture employed for estimating the range and Doppler of multiple targets and assume perfect timing synchronization at the RX side1. The exact spatial localization of the detected targets would require the adoption of multiple-input multiple output (MIMO) radar; however, this issue is out of the scope of this manuscript. 1The presence of a frequency offset does note represent a technical problem in this case, thanks to the fact that the TX and RX sides are co-located.
3 SEFDM-based TX/RX ISAC node RX comm. user Targets of interest Other scatterers Figure 1. Representation of the sensing scenario in a SEFDM-based ISAC system. In this case, sensing is accomplished by a wireless node that can establish a communication link with one or multiple users. The presence of targets of interests and other passive targets (scatterers) is explicitly indicated. Symbol mapping FrIDFT CP insertion Pulse shaping Sampling and CP removal Matched filtering DFT Symbol division Transmitter Data in Further processing Co-located receiver Point Targets Channel estimation Figure 2. Baseband architecture of the considered SEFDM-based co-located radar system. In the following, we also assume that: 1) The mth transmitted SEFDM symbol (with m= 0,1, ..., M −1) conveys the N-dimensional vector c(m) N≜ [c(m) 0, c(m) 1, ..., c(m) N−1]T, containing Nuchannel symbols (which belong to an Mc-ary constellation) and (N−Nu)zeros (associated with the suppressed subcarriers2). 2) Each SEFDM symbol contains a CP having size Ncp, lasting Tcp s, and whose presence guarantees the absence of ISI among adjacent SEFDM symbols. 3) The same subcarrier compression parameter 0< β ≤1, which is used to characterize the compression of the spectrum occupied by the Nsubcarriers compared to OFDM, is adopted for all the SEFDM symbols. 4) The spectrum P(f)of the pulse shaping filter, having impulse response (IR) p(t), used in the signal generation stage corresponds to a root of a raised cosine (RRC) filter with roll-off factor α. 5) At the RX side, a filter matched to p(t)is employed. The baseband architecture of the SEFDM-based radar system considered in this manuscript is illustrated in Fig. 2. Our derivation of the SEFDM signal model parallels that provided in [17, Sec. II-A] for a single OFDM symbol. The main difference is represented by the fact that the vector c(m) Nundergoes an order Nfractional inverse discrete Fourier transform (FrIDFT) of positive parameter β≤1(note that, if β= 1, SEFDM coincides with OFDM); this produces the N-dimensional vector x(m) N≜x(m) 0, x(m) 1, ..., x(m) N−1T=FH N,β c(m) N, (1) where FN,β is a square matrix of order N, whose (p, q)th element is equal to exp(−j2πpqβ/N)/√N. Then, x(m) Nis cyclically extended through a CP; for any m, the first Ncp elements of the resulting (N+Ncp)-dimensional vector ˘ x(m) Nare generated as ˘x(m) k=x(m) modN[k], with k∈ {−Ncp,−Ncp + 1, ..., −1}. 2Some subcarriers are deliberately suppressed to avoid out-of-band emissions, reduce interference with adjacent channels, and comply with pulse shaping constraints. Their use, as indicated in [17, Sec. II.A], notably simplifies receiver structure.
4 Under the above assumptions, the complex envelope of the transmitted signal for the considered SEFDM frame can be expressed as (e.g., see [10, Eq. (1)] for OFDM) ˜s(t) = M−1 X m=0 st−mT, x(m) N, (2) where (see [8, Eq. (2)] and the related comments) s(t, xN)≜ N−1 X k=−Ncp xkp(t−kTs/β), (3) T≜NTs/β and Tsis the channel symbol interval. The signal in (2) is transmitted over a multipath fading channel, where multiple propagation paths arise due to reflections from Lpoint-like objects in the environment. Among these reflectors, some may correspond to actual targets of interest, while others act as passive scatterers that contribute to the overall channel response (e.g., see Fig. 1). The resulting channel impulse response (CIR) is h(t, τ)≜ L−1 X l=0 ˜ hl(t, τ), (4) where ˜ hl(t, τ)≜alexp(j2πνlt)δ(τ−τl)(5) represents the CIR component associated with the lth target (characterized by the gain al, the delay τland the Doppler shift νl). In the following, we assume that: a) the CIR components are organized according to increasing delays, so that τ0and τL−1 represent the minimum and maximum delays, respectively; b) the parameters al,τland νldo not change over the entire frame (quasi-static channel). Note that, in this context, the delay τland Doppler shift νlcan be related to the physical parameters of the lth target, namely its range Rland radial velocity vl, since τl≜2Rl/c and νl≜2fcvl/c, where cdenotes the speed of light and fcthe carrier frequency. Let us derive now the expression of the matched filter output at the RX side when this filter is fed by r(t), i.e., the channel response to ˜s(t)(2). To simplify our developments, we first evaluate the response ˜rl(t)≜˜s(t)∗˜ hl(t, τ)∗p∗(−t), (6) which is obtained when the lth target only is active (so that the CIR is given by ˜ hl(t, τ)(5)) and channel noise is negligible; then, we account for the presence of multiple targets by summing over land adding the noise contribution. Following a similar approach as that illustrated in [16, Sec. II, Eqs. (4)-(7)], it can be proved that ˜rl(t) = M−1 X m=0 ¯rl,m(t)(7) for t∈(mNTs,(m+ 1)NTs), with m= 0,1, ..., M −1; here, ¯rl,m(t) = β √NTs ˜al N−1 X n=0 ¯x(m) nP(ϕn)P∗(ϕn−νl) expj2πϕn(t−τl)exp(j2πνlt), (8) represents the contribution of the mth SEFDM symbol, ϕn≜nβ∆f,∆f≜1/(NTs)is the subcarrier spacing,˜al≜ alexp(−j2πνlτl)and ¯x(m) nis the nth element of the order NDFT of x(m) N ¯ x(m) N≜DFTN[x(m) N] = ΞNx(m) N. (9) The signal ˜rl(t)(7) is sampled at the instant tm,˜n≜τL−1+ ˜nTs/β +mNT′, with ˜n= 0,1, ..., N +Ncp −1and T′≜ (Ts+Tcp)/β. Then, the resulting sequence {˜rl,m,˜n≜˜rl(tm,˜n)}undergoes CP removal; this produces the new sequence {rl,m,˜n}, with rl,m,˜n= ˘al N−1 X n=0 ¯x(m) nexp(−j2πnfτl) expj2πn ˜n Nd˜n(fνl) expj2πmfνl. (10) In the last expression, ˘al≜alexp(−j2πNfνlfτl)/√N, (11) fτl≜βτl−τL−1 NTs (12)
5 and fνl≜νl β∆f (13) are the complex gain, the normalized delay and normalized Doppler frequency, respectively, associated with the lth target, and d˜n(fνl)≜expj2π˜nfνl N(14) represents a phase rotation, proportional to the sample index ˜nand to the normalized Doppler frequency, associated with the lth target. For a given l, the samples {rl,m,˜n}are stored in the M×Nmatrix3 Rl≜[rl,m,˜n]=(X⊙Hl)ΞH NDl, (15) where X≜¯ x(0) N,¯ x(1) N, ..., ¯ x(M−1) NT(16) is the M×Nmatrix collecting all the channel symbols transmitted within a frame (see (9)), Hl≜[Hl,m,˜n]is an M×N channel matrix, with Hl,m,n ≜˘alexp(−j2πnfτl) expj2πmfνl(17) for any mand n,Dl≜diag(dl)is an N×Nmatrix and dl≜[d0(fνl), d1(fνl), ..., dN−1(fνl)]T(18) isaN-dimensional vector, whose ˜nth element is expressed by (14). Finally, the overall RX signal matrix is obtained by summing the contributions that originate from all the point targets and including the contribution of channel noise; this produces the M×Nmatrix R≜ L−1 X l=0 Rl+W= L−1 X l=0 (X⊙Hl)ΞH NDl+W, (19) where W≜[w[m, n]] is the M×Nmatrix representing the contribution of channel noise; in the following, it is assumed that the elements of Ware independent and identically distributed (i.i.d.) complex Gaussian random variables. B. Channel estimation Given the received signal model in (19), the radar channel, containing information about all reflectors in the considered scene, can be estimated at the ISAC receiver, which is assumed to be co-located with the transmitter (see Fig. 1). First of all, Rundergoes an order NDFT; this produces Y=R ΞN= L−1 X l=0 (X⊙Hl)ΞH NDlΞN+W ΞN. (20) The last result can be simplified if we assume that all the Doppler frequencies {νl}do not exceed the SEFDM subcarrier spacing, i.e., that |νl|< β∆ffor any l; this is equivalent to assuming a limited channel variations, as usually done in the study of OFDM-based sensing. In fact, under the last assumption, ΞH NDlΞN∼ =INfor any l, so that (20) can be rewritten as Y∼ =X⊙H+¯ W, (21) where H≜ L−1 X l=0 Hl(22) represents the overall channel matrix and ¯ W≜WΞN. Since the receiver is assumed to be co-located with the transmitter, the data symbol matrix X(16) is known; consequently, an estimate of H(22) can be obtained from (21) by evaluating4 ˆ H=Y⊘X∼ =H+˘ W, (23) where ˘ W≜¯ W⊘Xis an M×Nnoise matrix. As it can be easily inferred from (17) and (22), the elements of ˆ H(23) form a 2D sequence consisting of Ldistinct complex tones superimposed with noise. The parameters of each tone provide information about the range, the radial speed and the radar cross section of a specific point target. Consequently, radar sensing can be performed by applying an algorithm for 2D harmonic retrieval to the matrix ˆ H(23), that collects the available set of noisy measurements. 3The dependence of the matrix Hlon the target parameters (˘al, fτl, fνl)is not explicitly shown in (15) and in the following to ease notation. 4A conceptually similar symbol division approach has been proposed in the context of OFDM-based radar systems (e.g., see [15, Eq. (5)]).
6 It is worth noting that, after the DFT operation, in (20), the transformed noise matrix ¯ Wstill contains i.i.d. Gaussian noise samples due to the unitary nature of the DFT operator. However, when computing the element-wise division in (23), the entries of ¯ Ware divided by the known data symbols in X, which belong to a modulation constellation (e.g., phase-shift keying, PSK, or quadrature amplitude modulation, QAM). As a consequence, the noise samples in the resulting matrix ˘ Ware no longer i.i.d. and their variances depend on the modulus of the corresponding entries in X. On the one hand, when a QAM modulation is used, the amplitude of the symbols varies, leading to a non-uniform noise power across the matrix; this is exacerbated in the case of SEFDM due to the intentional ICI. This implies that the classical ML estimation metric, which typically assumes i.i.d. noise, becomes suboptimal in both OFDM and SEFDM. On the other hand, when a constant-modulus constellation such as PSK is employed, all the elements of Xhave the same magnitude, and the samples of ˘ Wretain a uniform variance, making the AWGN assumption valid and preserving the optimality of the classical ML metric when β= 1 (i.e., in the OFDM case) only. A novel technique for accomplishing this last task is described in the following section, under the assumption of i.i.d. noise samples affecting the channel estimates (23). It is worth pointing out that the adoption of a CP in SEFDM, similarly to OFDM, has an important impact on both SE and estimation accuracy. In fact, on the one hand, the CP reduces the overall SE. On the other hand, however, it makes the adoption of the Fourier-based representation (8) possible and eliminates ISI, thus allowing the straightforward channel estimation procedure5expressed by (23). Conversely, if the CP was not employed, (8) would not hold and adjacent blocks would interfere, so that the simple estimation method in (23) would no longer be applicable; therefore, substantially more complicated estimation methods would be required. III. APPROXIMATE MAXIMUM LIKELIHOOD ESTIMATION OF CHANNEL PARAMETERS This section is organized into two parts. In the former, we develop a novel algorithm for estimating the set SH≜ {(˘al, fνl, fτl); l= 0,1, ..., L −1}, that collects the parameters of the L2D complex exponentials contributing to the channel matrix H(22). In the second part, we analyze the computational complexity of the proposed method. A. Derivation of the proposed algorithm To begin, we rewrite the model (23) in vector form as ˆ h≜vec( ˆ H)=h+˘ w, (24) where h≜vec(H)and ˘ w≜vec( ˘ W)are (MN)-dimensional column vectors. Based on (17) and (22), the vector hcan be expressed as h(a,fτ,fν)=B(fτ,fν)a; (25) here, B(fτ,fν)≜C∗ N(fτ)⊛CM(fν)(26) is an (MN)×Lmatrix, Cp(f)≜[¯ b0(f),¯ b1(f), ..., ¯ bp−1(f)]T(27) isap×Lmatrix for any positive integer p, ¯ bx(f)≜[¯ bx(f0),¯ bx(f1), ...,¯ bx(fL−1)]T(28) is an L-dimensional column vector and ¯ bx(fl)≜exp(j2πxfl)for any l. Moreover, a≜[˘a0,˘a1, ..., ˘aL−1]T, (29) fτ≜[fτ0, fτ1, ..., fτL−1]T(30) and fν≜[fν0, fν1, ..., fνL−1]T(31) are the L-dimensional column vectors collecting the complex gains, the normalized delays and normalized Doppler frequencies, respectively, that characterize the considered channel. Given the estimate ˆ h(24) of h, our goal is to determine its inner structure by estimating Land SH; this is equivalent to solving the problem of target detection and estimation in a SEFDM-based radar system, in which the lth target is modeled as a point target with parameters (˘al, fτl, fνl). For this reason, if Lis known, the problem of target estimation can be formulated as a ML estimation problem. Specifically, given the observation vector ˆ h, we aim to estimate the unknown target parameters (a,fτ,fν), under the assumption that the noise vector ˘ win (24) is a realization of a zero-mean circularly symmetric complex Gaussian random vector with independent 5Note this procedure makes an accurate recovery of the channel matrix H(22) possible in the absence of channel noise.
7 entries, having zero mean and variance σ2 w. Given this assumption, the conditional probability density function (PDF) of ˆ h, conditioned on the trial values (˜ a,˜ fτ,˜ fν)of the mentioned unknown parameters, is p(ˆ h|˜ a,˜ fτ,˜ fν) = 1 (πσ2 w)MN exp−1 σ2 w ˆ h−h(˜ a,˜ fτ,˜ fν) 2. (32) The ML estimate (ˆ a,ˆ fτ,ˆ fν)of (a,fτ,fν)is therefore achieved by maximizing the likelihood function (32), as (ˆ a,ˆ fτ,ˆ fν)≜arg max ˜ a,˜ fτ,˜ fν p(ˆ h,|,˜ a,˜ fτ,˜ fν). (33) Equivalently, since the exponential function, in (32), is monotonically increasing, the ML estimation problem can be reformulated as the minimization of the squared error appearing in the same function. This easily leads to (ˆ a,ˆ fτ,ˆ fν)≜arg min ˜ a,˜ fτ,˜ fνL(ˆ h|˜ a,˜ fτ,˜ fν), (34) where L(ˆ h|˜ a,˜ fτ,˜ fν)≜Trnˆ h−h(˜ a,˜ fτ,˜ fν)ˆ h−h(˜ a,˜ fτ,˜ fν)Ho(35) is a log-likelihood (LL) cost function for the optimization problem (34). Given ˜ fτand ˜ fν, the minimum of L(ˆ h|˜ a,˜ fτ,˜ fν)with respect to ˜ ais obtained by 1) substituting the expression of h(·,·,·)(25) in the RHS of (35), 2) taking the derivative of the resulting expression with respect to ˜ a, 3) setting it to zero and 4) solving for ˆ a=˜ a; this results in ˆ a=B†(˜ fτ,˜ fν)ˆ h. (36) Unfortunately, given ˜ a=ˆ a, the minimization of (35) with respect to ˜ fτand ˜ fνdoes not lead to a closed-form solution. However, if a coarse estimate of both ˜ fτand ˜ fνis available, a Newton-based method can be adopted to solve the optimization problem in (34). This approach requires evaluating the gradient vector and the Hessian matrix of L(35) with respect to the vector ˜ f≜[˜ fT τ˜ fT ν]T. Proposition 1: The gradient vector and the Hessian matrix of Lcan be expressed as6 ∇L(˜ a,˜ f) = −2hℜ˜ a∗⊙˙ BH ˜ fτ(ˆ h−B˜ a)T,ℜ˜ a∗⊙˙ BH ˜ fν(ˆ h−B˜ a)TiT , (37) and as ¨ HL(˜ a,˜ f)≜T˜ fτ,˜ fτU˜ fτ,˜ fν U˜ fτ,˜ fνV˜ fν,˜ fν, (38) respectively; here7, T˜ fτ,˜ fτ≜"∂2L ∂˜ fτl∂˜ fτl′#= 2ℜn(˜ a˜ aH)⊙˙ BH ˜ fτ ˙ B∗ ˜ fτ−˜ a∗ˆ h−B˜ aH¨ B˜ fτ,˜ fτo, (39) U˜ fτ,˜ fν≜"∂2L ∂˜ fτl∂˜ fνl′#= 2ℜn(˜ a˜ aH)⊙˙ BH ˜ fτ ˙ B∗ ˜ fν−˜ a∗ˆ h−B˜ aH¨ B˜ fτ,˜ fνo (40) and V˜ fν,˜ fν≜"∂2L ∂˜ fνl∂˜ fνl′#= 2ℜn(˜ a˜ aH)⊙˙ BH ˜ fν ˙ B∗ ˜ fν−˜ a∗ˆ h−B˜ aH¨ B˜ fν,˜ fνo (41) are L×Lmatrices. Moreover, ˙ B˜ fxis an (MN)×Lmatrix, whose lth column contains the partial derivative of the lth column of B(see (26) with (fτ,fν)=(˜ fτ,˜ fν)), evaluated with respect to ˜ fx. Similarly, ¨ B˜ fx,˜ fyis a (MN)×Lmatrix, whose lth column contains the second-order derivative of the lth column of B(26) evaluated with respect to ˜ fxand ˜ fy. Note that, thanks to the linearity property of the derivative, the equality ¨ B˜ fx,˜ fy=¨ B˜ fy,˜ fxholds. Proof: The results illustrated in this proposition are proven in Appendix A . Given the previous results, an iterative procedure for estimating the parameters a,fτand fνcan be formulated as follows: 1) Initialization – The initial values ˆ a(0),ˆ f(0) τand ˆ f(0) νare obtained through the 2D periodogram method8, that also generates an estimate ˆ Lof L. This method is based on a peak search on the 2D spectrum of ˆ H(23) and produces ˆ Lcouples of 2D frequencies, associated with the delay and Doppler of ˆ Ltargets. Here, we assume that all point targets are sufficiently spaced along the delay and Doppler dimensions, so that the peaks appearing in the 2D spectrum evaluated for the received signal 6The dependence of the matrix B(26) on the variables (˜ fτ,˜ fν)is not explicitly shown in the following to ease reading. 7Note that the dependence of the matrices T˜ fτ,˜ fτ,U˜ fτ,˜ fνand V˜ fν,˜ fνon the variable ˜ ais not explicitly shown to ease notation. 8The 2D periodogram method is dubbed 2D-FFT in the following.
8 can be easily identified without ambiguity. Under this condition9, an estimate10 ˆ Lof Lcan be obtained by adopting a proper detection threshold ϵinit in the evaluation of the local maxima of the above mentioned 2D spectrum. Note that this threshold is typically selected in a way to minimize the detection of false (i.e., ghost) spectral components; for instance, the generalized likelihood ratio test, GLRT, can be been employed [18, Par. 4.6.1]). Then, the ˆ Lestimated 2D frequencies are organized in the (2ˆ L)-dimensional vector ˆ f(0) = [(ˆ f(0) τ)T(ˆ f(0) ν)T]T. Finally, the iteration index iis set to 1. 2) Refinement procedure – An iterative procedure, consisting of the following three steps, is executed. 2a) Complex amplitude update – The new complex amplitude vector ˆ a(i)is computed through (36) with (˜ fτ,˜ fν) = (ˆ f(i−1) τ,ˆ f(i−1) ν), where ˆ f(i−1) τand ˆ f(i−1) νdenote the estimates of fτand fν, respectively, available at the end of the (i−1)th iteration. 2b) Frequency update – The new estimate ˆ f(i)=ˆ f(i−1) −µ¨ H−1 Lˆ a(i),ˆ f(i−1)∇Lˆ a(i),ˆ f(i−1), (42) is computed; here, ˆ f(i)= [(ˆ f(i) τ)T(ˆ f(i) ν)T]Tand µis a real positive parameter representing the step size of the proposed method. In our simulations, the step size µhas been set to one11. 2c) Convergence test – The quantity ∆L(i)≜Lˆ h|ˆ a(i),ˆ f(i) τ,ˆ f(i) ν−Lˆ h|ˆ a(i−1),ˆ f(i−1) τ,ˆ f(i−1) ν(43) is evaluated (see (35)). Then, this last value is compared with a proper threshold12 ϵit.If∆L(i)> ϵit, the iteration index iis increased by one and, if i<Nit, steps 2a)-2b)-2c) are repeated. Otherwise, the algorithm stops, and the quantities ˆ a(i),ˆ f(i) τ and ˆ f(i) νare taken as output. Algorithm 1: Newton-based multiple cisoid refiner (NMCR) Input: The vector ˆ h(24), the detection threshold ϵinit, the oversampling factors (LD, Lr)for the 2D-FFT method, the parameter ϵit, the step size µand the overall number of iterations Nit. 1Initialization: Evaluate the initial estimates ˆ a(0),ˆ f(0) τ,ˆ f(0) νand ˆ Lby applying the 2D-FFT method to ˆ h(24) . Then, set the iteration index ito 1. 2Refinement: 2aComplex amplitude update: compute ˆ a(i)through (36) with (˜ fτ,˜ fν)=(ˆ f(i−1) τ,ˆ f(i−1) ν). 2bFrequency update: Compute ˆ f(i)= [(ˆ f(i) τ)T(ˆ f(i) ν)T]Tthrough (42). 2cConvergence test: Evaluate the quantity ∆L(i)(43); then, if ∆L(i)> ϵit, the iteration index iis increased by one and, if i < Nit, steps 2a)-2b)-2c) are repeated. Output: The estimates ˆ a(i),ˆ f(i) τ,ˆ f(i) νand ˆ Lof a,fτ,fνand L, respectively. The NMCR algorithm is summarized in Algorithm 1. B. Computational complexity analysis It is not difficult to show that the overall computational cost of a single iteration of the refinement procedure of the NMCR algorithm in the presence of Ltargets can be expressed as O((L+L2)M N). The proposed procedure requires an initial estimate of the frequencies fτand fν. If the 2D-FFT method is used to compute this estimate, the complexity of the initialization is O(N2D-FFT), with (e.g., see [16, Sec. III-A.1, Eq. (26)]) N2D-FFT =M0N0log2(M0N0)+M0N0+Llog2(L); (44) here, M0≜MLD,N0≜NLr, whereas LDand Lrare the so-called oversampling factors that are selected for the discrete symplectic Fourier transform (DSFT) operation13 executed along the rows and columns of ˆ H(23), respectively. Therefore, the overall computational cost of the NMCR algorithm is O(NNMCR), with NNMCR ≜N2D-FFT +Nit(L+L2)M N. (45) 9If this condition is not met, closely spaced targets may cause partially overlapped peaks, potentially resulting in an incorrect estimate of ˆ Lor a degraded initialization accuracy. 10For notational simplicity, in the following we assume ˆ L=Lwithout loss of generality. It should be noted, however, that the proposed algorithm operates on the available estimate ˆ L, so its applicability is preserved even when ˆ L=L. 11In principle, one could reduce the value of the step size µ(for instance, it could be halved) if no significant decrease in the value of L(35) is observed over two consecutive iterations. 12The threshold ϵit can be selected using the same approach adopted for the threshold employed for the initialization step, i.e., ϵinit. 13Further details about the 2D-FFT method and its inherent DSFT processing can be found in [16, Sec. III-A.1].
9 IV. NUMERICAL RESULTS In this section, we assess the accuracy of the NMCR algorithm in six different scenarios and we compare it with that achieved by the following four algorithms: 1) the 2D-FFT method [9]; 2) a refinement procedure based on interpolation of the 2D periodogram [19, Sec. IV-A]; 3) the approximate ML method, based on alternating projections (APs), recently proposed in [14]; 4) the expectation maximization (EM) algorithm [14]. In the following, the last three methods are denoted 2D-FFTi, AP-ML and EM, respectively. Note that such methods have originally been developed in the context of OFDM, but can be easily adapted to SEFDM. The characteristics of the considered scenarios can be summarized as follows (the ith scenario is denoted Siin the following): S1) This is characterized by a single target (i.e., L= 1) having unitary amplitude A0, and whose range R0and velocity v0 are random variables uniformly distributed in [Rmin, Rmax]and [vmin, vmax], respectively; in our simulations, [Rmin, Rmax] = [5,45] m and [vmin, vmax] = [5,50] m/s have been selected. S2) This is characterized by L= 2 targets. The parameters (A0, R0, v0)of the first target are generated in the same way as those given for S1, whereas those of the second target, namely (A1, R1, v1)have been selected in a way that A1= 0.8+0.2 ¯α A0, R1=R0+¯ R Rbin and v1=v0+ ¯v vbin, respectively; here, ¯αis a random variable uniformly distributed over the interval [0,1],¯ Rand ¯vrepresent the normalized spacing in range and Doppler, respectively (both have been set to 3in our simulations referring to S2), and Rbin ≜c 2Nβ∆f ,vbin ≜c β 2MfcTs (46) denote the range and velocity resolutions, respectively, for the considered SEFDM-based radar system. S3) This is characterized by L= 5 targets. The parameters (A0, R0, v0)of the first target have been generated in the same way as S1; however, [Rmin, Rmax] = [10,35] m and [vmin, vmax] = [30,50] m/s have been chosen in the generation of R0and v0. Moreover, the following choices have been made for the lth target (with l= 1,2,3,4): 1) Rl=R0+ 2l(¯α−0.5)Rbin; 2) vl=v0; 3) Al, being generated according to Swerling-3model (e.g., see [20, Eq. (7)]), follows a chi-square distribution with a mean value equal to A0. S4) This is characterized by L= 3 targets. The parameters of the first target have been generated in the same way as S1. Moreover, the following choices have been made for the lth target (with l= 1,2): 1) Al=A0= 1; 2) the target ranges and velocities have been generated in the same way as S2, but ¯ R= 1.25 and ¯v= 1.75 have been selected. S5) This is characterized by a variable number of targets L∈ {1,2, ..., 8}. The parameters of the first target have been generated in the same way as S3, whereas those of the lth target (with l= 1,2, ..., L −1) have been produced in the same way as S2, but ¯ R= ¯v= 1 have been chosen. S6) This scenario is characterized by a couple of targets (i.e., L= 2), whose physical parameters are generated in the same way as S2; however, we have that: a) A1= 0.75, b) ¯ R= 2 and c) ¯v= 2. Note that the first three scenarios are characterized by a variable SNR ∈[−20,10] dB14, whereas the SNR has been set to −5dB in both S4 and S5, and to 0dB in S6. The selection of the scenarios described above can be motivated as follows: 1) In S1, which is characterized by L= 1, target detection and estimation are not affected by the potential spectral leakage originating from the presence of multiple targets. 2) In S2, multiple targets are present, but leakage is negligible, the targets being adequately spaced. 3) S3 is useful to assess the estimation accuracy in the presence of an extended scatterer. In fact, in this case, the targets have a common centroid (characterized by the parameters (A0, R0, v0)), share the same velocity and are closely spaced in range. Moreover, the target amplitudes follow a model typically used for extended targets. 4) S4 and S5 are characterized by multiple closely spaced targets, each with a fixed SNR. In our analysis, S4 has been specifically designed to examine the convergence behavior of the proposed NMCR algorithm and compare it with other iterative techniques, namely the AP-ML and EM algorithms. To this end, the dependence of the estimation accuracy on the number of iterations, employed in the frequency refinement step of each method, has been assessed. Therefore, readers should keep in mind that S4 provides insight into the stability and robustness of the considered iterative procedures, highlighting how quickly and reliably each algorithm approaches its steady-state performance. In contrast, S5 has been adopted to assess the dependence of the computational complexity of each method on the number of targets. 5) S6 has been considered to assess the impact of the compression factor βon sensing performance. A change in this parameter modifies the normalized delay and Doppler frequency of each target, but does not alter the normalized frequency resolutions, which are only determined by the parameters Mand N. In our simulations, we have also considered an additional scenario (denoted as S0), whose exclusive scope is the analysis of the detection capability of the GLRT method employed in the initialization step of all the considered algorithms under different SNR conditions, number of targets, and target velocities. In fact, for this scenario, we have assessed the detection probability PD, i.e., the probability15 that the estimated number of targets ˆ Lcoincides with the true one L. The following assumptions have been made about the targets: 1) L∈ {1,3,5}; 2) the parameters of the first target are generated in the same 14In all the considered scenarios, the impact of the CP on the evaluation of the effective SNR has been considered. 15The GLRT threshold ϵinit has been selected in a way to ensure a constant false alarm probability Pfa = 10−3.
16 Possible directions for future research include: 1) the extension of the proposed framework to MIMO radar architectures, enabling joint range–Doppler–angle estimation and full target localization; 2) hybrid waveform design for ISAC; 3) the adoption of alternative non-orthogonal modulation formats that are more robust to doubly-selective channels. APPENDIX A. Gradient and Hessian of the LL function In this appendix, the derivation of the gradient (37) and of the Hessian matrix (38) for the cost function L(35) is sketched. Before doing so, some basic properties of the Dirac vectors as well as of the stacking operations are illustrated. To this aim, we define: a) the L-dimensional vectors x≜[x0, x1, ..., xL−1]T,y≜[y0, y1, ..., yL−1]Tand z≜[z0, z1, ..., zL−1]T(the lth element is denoted xl,yland zl, respectively); b) the L×Lmatrix A≜[al,l′](its (l, l′)th element is denoted al,l′); c) the L-dimensional vector δL,l, called Dirac vector, having its lth element equal to unity and its remaining (L−1) elements equal to zero. Given these definitions, the following properties can be proved (here, Piindicates the ith property): •P1:xTδL,l =xl •P2:δL,l AδL,l′=al,l′ •P3:ifzl=xlylfor any l, then z=x⊙y •P4:ifal,l′=xlyl′, for any land l′, then A=xyT Our derivations start with substituting the expression of h(·,·,·)(25) in the RHS of (35) and evaluating the partial derivative of the resulting expression with respect to the scalar frequency ˜ Fτl(with l= 0,1, ..., L −1); this yields ∂L ∂˜ fτl =−2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fτˆ h−B˜ ao, (47) where ˙ BH ˜ fτis an (MN)×Lmatrix, whose lth column contains the partial derivative of the lth column of B(see (26) with (fτ,fν) = (˜ fτ,˜ fν)), with respect to ˜ fτl. The partial derivative of L(35) with respect to the scalar frequency ˜ fνl(with l= 0,1, ..., L −1) can be computed in a similar way; this results in ∂L ∂˜ fνl =−2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fνˆ h−B˜ ao, (48) where ˙ B˜ fνis an (MN)×Lmatrix, whose lth column contains the partial derivative of the lth column of B(see (26) with (fτ,fν)=(˜ fτ,˜ fν)), with respect to ˜ fνl. Given the partial derivatives in (47) and (48), the gradient vector ∇L(˜ a,˜ f)(having length equal to 2L) can be obtained, since ∇L(˜ a,˜ f)≜"∂L ∂˜ fτ0 ,∂L ∂˜ fτ1 , ..., ∂L ∂˜ fτL−1 ,∂L ∂˜ fν0 ,∂L ∂˜ fν1 , ..., ∂L ∂˜ fνL−1#T . (49) Exploiting the property P1 in both (47) and (48) and the property P3 in (49) easily leads to (37). The evaluation of the (2L)×(2L)Hessian matrix ¨ HL(˜ a,˜ f)requires computing the partial derivatives of L(35) with respect to the couple of variables (˜ fτl,˜ fτl′),(˜ fτl,˜ fνl′),(˜ fνl,˜ fτl′)and (˜ fνl,˜ fνl′). Following the same line of reasoning as the one illustrated for the evaluation of the gradient of L, it can be easily shown that ∂2L ∂˜ fτl∂˜ fτl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fτ,˜ fτ(ˆ h−B˜ a)o+ 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fτ ˙ B˜ fτ δL,l′δT L,l′˜ ao, (50) ∂2L ∂˜ fτl∂˜ fνl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fτ,˜ fν(ˆ h−B˜ a)o+ 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fτ ˙ B˜ fν δL,l′δT L,l′˜ ao, (51) ∂2L ∂˜ fνl∂˜ fτl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fν,˜ fτ(ˆ h−B˜ a)o+ 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fν ˙ B˜ fτ δL,l′δT L,l′˜ ao(52) and ∂2L ∂˜ fνl∂˜ fνl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fν,˜ fν(ˆ h−B˜ a)o+ 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fν ˙ B˜ fν δL,l′δT L,l′˜ ao. (53) Given equations (50)-(53), the (2L)×(2L)Hessian matrix of ¨ HL(˜ a,˜ f)can be easily put in the form (38), by exploiting: 1) P1 for all the terms ˜ aHδL,l and δT L,l′˜ ain (50)-(53); 2) P1 for the terms δT L,l′¨ BH ˜ fτ,˜ fτ y,δT L,l′¨ BH ˜ fτ,˜ fν y,δT L,l′¨ BH ˜ fν,˜ fτ yand δT L,l′¨ BH ˜ fν,˜ fν y, in (50), (51), (52) and (53), respectively (here, y≜ˆ h−B˜ a); 3) P2 for the terms δT L,l ˙ BH ˜ fτ ˙ B˜ fτ δL,l′,δT L,l ˙ BH ˜ fτ ˙ B˜ fν δL,l′,δT L,l ˙ BH ˜ fν ˙ B˜ fτ δL,l′and δT L,l ˙ BH ˜ fν ˙ B˜ fν δL,l′, in (50), (51), (52) and (53), respectively. 4) P3 and P4 in the resulting expressions.
17 B. Derivation of the Cramér-Rao Lower Bounds In this appendix, the CRLB for the joint estimation of the complex gains a(29), the normalized delays fτ(30) and the normalized Doppler frequencies fν(31) associated with Ltargets is derived. Note that the CRLBs referring to the considered estimation problem have already been derived for OFDM-based radar systems (e.g., see [12]); however, as far we know, the same derivation for SEFDM-based co-located radar systems has never appeared so far. First of all, let us consider the model (24), that refers to Ldistinct targets, and define the trial vectors ˜ a,˜ fτand ˜ fνin a similar way as a(29), fτ(30) and fν(31), respectively. The CRLBs we are interested in refer to the ML estimation problem (34). If we assume that the elements of the noise vector ˘ w, in (24), are Gaussian, mutually independent and have zero mean and variance σ2 w, the CRLBs of all the parameters of interest are represented by the diagonal elements of the matrix V=σ2 wΦ−1, (54) where Φ≜[Φi,j] = 2ℜ ∂˜ h ∂˜ f ∂˜ h ∂˜ f!H (55) is the (2L)×(2L)Fisher information matrix (FIM) computed for the (MN)-dimensional row vector ˜ h≜hT(˜ a,˜ fτ,˜ fν)(see (25)). Note that, in (55), ˜ f≜[˜ fT τ,˜ fT ν]Tisa(2L)-dimensional column vector. Moreover, it is not difficult to prove that ∂˜ h ∂˜ f= [P,Q]T, (56) where P≜[˜ P0,˜ P1, ..., ˜ PL−1]and Q≜[˜ Q0,˜ Q1, ..., ˜ QL−1]are (MN)×Lmatrices containing the partial derivatives with respect to the delay and Doppler parameters {fτl}and {fνl}, respectively. The lth column of Pand Qcan be expressed as (see (28)) ˜ Pl=−˘alΥN⊙¯ b∗ N(fτl)⊛¯ bM(fνl), (57) ˜ Ql= ˘al¯ b∗ N(fτl)⊛ΥM⊙¯ bM(fνl), (58) respectively; finally, ΥX≜[0, j2π, ..., j2π(X−1)]T, for any integer X. It is worth emphasizing that the FIM in (55) includes cross-terms requiring the evaluation of partial derivatives with respect to heterogeneous parameters. In fact, substituting the RHS of (56) in that of (55) yields Φ= 2ℜPHP PHQ QHP QHQ (59) The presence of the above-mentioned cross-terms results in a dependence of the last matrix, and consequently, of the matrix V(54), on the relative spacing between the normalized delays and Doppler frequencies of the Ltargets. Consequently, even when the actual range and Doppler values of the targets are kept fixed, any change in the compression factor β(which scales the normalized frequency grid) affects the values of the normalized frequencies fτand fν, and thus modifies the structure of the FIM and the resulting CRLBs. This justifies the β-dependent behavior of the considered estimation bounds, as discussed in the main text. In our work, all the above-mentioned CRLBs have been evaluated numerically on the basis of (54)-(59); in doing so, the expected value of each random variable in the set {(˘al, fτl, fνl); l= 0,1, ..., L −1}has been used. REFERENCES [1] J. A. Zhang et al., “An Overview of Signal Processing Techniques for Joint Communication and Radar Sensing,” IEEE J. Sel. Topics Signal Process., vol. 15, no. 6, pp. 1295–1315, Nov. 2021. [2] Z. Wei et al., “Integrated Sensing and Communication Signals Toward 5G-A and 6G: A Survey,” IEEE Internet Things J., vol. 10, no. 13, pp. 11 068– 11 092, Jul. 2023. [3] I. Darwazeh and M. Rodrigues, “A Spectrally Efficient Frequency Division Multiplexing Based Communications System,” in Proceedings of the 8th International OFDM-Workshop (InOWo’03), Sep. 2003. [4] S. Isam, I. Kanaras, and I. Darwazeh, “A Truncated SVD approach for fixed complexity spectrally efficient FDM receivers,” in 2011 IEEE Wireless Commun. and Netw. Conf. (WCNC), Mar. 2011, pp. 1584–1589. [5] S. Isam and I. Darwazeh, “Design and Performance Assessment of Fixed Complexity Spectrally Efficient FDM Receivers,” in 2011 IEEE 73rd Veh. Tech. Conf. (VTC Spring), May 2011, pp. 1–5. [6] B. Yu et al., “Channel equalisation and data detection for SEFDM over frequency selective fading channels,” IET Communications, vol. 12, no. 18, pp. 2315–2323, Oct. 2018. [7] Y. Ma et al., “A Low-Complexity Receiver for Multicarrier Faster-Than-Nyquist Signaling Over Frequency Selective Channels,” IEEE Commun. Lett., vol. 24, no. 1, pp. 81–85, Jan. 2020. [8] M. Mirabella, P. Di Viesti, and G. M. Vitetta, “A Novel Message Passing Algorithm for Soft-Output Detection in Faster-than-Nyquist Multicarrier Systems,” in 2025 IEEE 26th Int. Workshop on Signal Process. and Artificial Intell. for Wireless Commun. (SPAWC), 2025, pp. 1–5. [9] C. Sturm and W. Wiesbeck, “Waveform Design and Signal Processing Aspects for Fusion of Wireless Communications and Radar Sensing,” Proc. IEEE, vol. 99, no. 7, pp. 1236–1259, Jul. 2011.
18 [10] S. Mercier, S. Bidon, D. Roque, and C. Enderli, “Comparison of Correlation-Based OFDM Radar Receivers,” IEEE Trans. Aerosp. Electron. Syst., vol. 56, no. 6, pp. 4796–4813, Dec. 2020. [11] L. Zheng and X. Wang, “Super-Resolution Delay-Doppler Estimation for OFDM Passive Radar,” IEEE Trans. Signal Process., vol. 65, no. 9, pp. 2197–2210, May 2017. [12] R. Xie, D. Hu, K. Luo, and T. Jiang, “Performance Analysis of Joint Range-Velocity Estimator With 2D-MUSIC in OFDM Radar,” IEEE Trans. Signal Process., vol. 69, pp. 4787–4800, Aug. 2021. [13] U. Singh, R. Mitra, V. Bhatia, and A. Mishra, “Target range estimation in OFDM radar system via kernel least mean square technique,” in Int. Conf. on Radar Sys. (Radar 2017). IET, 2017, p. 44. [14] F. Zhang, Z. Zhang, W. Yu, and T.-K. Truong, “Joint Range and Velocity Estimation With Intrapulse and Intersubcarrier Doppler Effects for OFDM-Based RadCom Systems,” IEEE Trans. Signal Process., vol. 68, pp. 662–675, Jan. 2020. [15] Y. L. Sit, C. Sturm, and T. Zwick, “Doppler estimation in an OFDM joint radar and communication system,” in 2011 German Microwave Conf., Mar. 2011, pp. 1–4. [16] M. Mirabella, P. Di Viesti, A. Davoli, and G. M. Vitetta, “Deterministic Signal Processing Techniques for OFDM-Based Radar Sensing: An Overview,” IEEE Access, vol. 11, pp. 68 872–68 889, Jul. 2023. [17] M. Mirabella, P. Di Viesti, and G. M. Vitetta, “On the Use of a Two-Dimensional Cyclic Prefix in OTFS Modulation and Its Implications,” IEEE Open J. Commun. Soc., vol. 5, pp. 3340–3367, 2024. [18] B. Ottersten, M. Viberg, P. Stoica, and A. Nehorai, “Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing,” in Radar array processing. Springer, 1993, pp. 99–151. [19] M. Mirabella, P. Di Viesti, and G. M. Vitetta, “Deterministic Algorithms for Four-Dimensional Imaging in Colocated MIMO OFDM-Based Radar Systems,” IEEE Open J. Commun. Soc., vol. 4, pp. 1516–1543, Jul. 2023. [20] D. Shnidman, “Expanded Swerling target models,” IEEE Trans. Aerosp. Electron. Syst., vol. 39, no. 3, pp. 1059–1069, Jul. 2003.