# On the Effectiveness of OTFS for Joint Radar and Communication

## Authors

• 3 publications
• 18 publications
• 118 publications
• 6 publications
02/04/2019

### Performance Analysis of Joint Radar and Communication using OFDM and OTFS

We consider a joint radar estimation and communication system using orth...
05/24/2021

### Frequency Permutations for Joint Radar and Communications

This paper presents a new joint radar and communication technique based ...
06/30/2020

### A Novel Bistatic Joint Radar-Communication System in Multi-path Environments

Radar detection and communication can be operated simultaneously in join...
09/17/2020

### Hybrid Digital-Analog Beamforming and MIMO Radar with OTFS Modulation

Motivated by future automotive applications, we study some joint radar t...
09/01/2021

### A Novel ISAC Transmission Framework based on Spatially-Spread Orthogonal Time Frequency Space Modulation

In this paper, we propose a novel integrated sensing and communication (...
07/30/2021

### Combined Radar and Communications with Phase-Modulated Frequency Permutations

This paper focuses on the combined radar and communications problem and ...
09/10/2019

### MAJoRCom: A Dual-Function Radar Communication System Using Index Modulation

Dual-function radar communication (DFRC) systems implement both sensing ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

A key-enabler of high-mobility networks is the ability of a node to continuously track its dynamically changing environment (state) and react accordingly. Although state sensing and communication have been designed separately in the past, power and spectral efficiency and hardware costs encourage the integration of these two functions, such that they are operated by sharing the same frequency band and hardware (see e.g. [1]). Motivated by emerging vehicular applications (V2X) [2], we consider a joint radar and communication system where a radar-equipped transmitter wishes to transmit information to a target receiver and simultaneously estimate the parameters of that receiver such as range and velocity.

Such a communication setup has been extensively studied in the literature (see [3, 4, 5, 6, 7, 8, 9, 10, 11]

and references therein). Existing works on joint radar and communications can be roughly classified into two classes. The first class considers some resource-sharing approach, such that time, frequency, or space resources are split into either radar or data communication (e.g., see

[7, Section III. A, C] and references therein). The second class uses a common waveform for both radar and communication. This approach includes information-embedded radar waveforms (e.g., see [7, Section III. D], [1, Section IV. A], [10], and references therein) as well as the direct usage of standard communication waveforms applied to radar detection (e.g., see [1, Section IV. B], [1, Section III. B], [3, 12, 13, 6, 14]). It is worth noticing that a synergistic waveform design for joint radar and communications yields a significant potential gain compared to the resource-sharing approach, as demonstrated in a simplified albeit representative information theoretic framework in [15]. It is also worthwhile to notice that typical V2X channels are characterized by a relatively small number of discrete multi-path components corresponding to line-of-sight (LoS) propagation, ground reflection, and some specular reflections on surrounding (e.g., metal) surfaces, each of which is characterized by its own, possibly large, Doppler frequency shift (e.g., see [16, 17] and references therein). Motivated by this fact, we focus on a joint radar and communication system using orthogonal time frequency space (OTFS) modulation (see [18, 19] and references therein), as this modulation format provides inherent robustness to Doppler shifts and is naturally suited to sparse channels in the delay-Doppler domain [20, 21].

The first part of this paper focuses on the suitability of OTFS for joint radar parameter estimation and data communications, by considering the mean square error (MSE) of range and velocity estimation and the achievable rate with Gaussian inputs for a simple point-to-point communication scenario. We extend our preliminary work [22], which considered only LoS propagation, to the case of a channel with multiple paths, one of which is the LoS. We propose an efficient approximated Maximum Likelihood (ML) algorithm to estimate the range and velocity of the target from the backscattered signal and derive the corresponding Cramér-Rao lower bound (CRLB). Our numerical examples inspired by the parameters of IEEE 802.11p demonstrate that digital multi-carrier modulation formats such as OTFS and orthogonal frequency division multiplexing (OFDM) yield as accurate estimation performance as frequency modulated continuous wave (FMCW), one of the typical automotive radar waveforms [4], while achieving significant communication rates. These results suggest that joint radar and communication can be effectively operated without compromising neither the achievable data rates nor the radar performance.

The second part of the paper addresses more specifically the soft-output symbol detection for OTFS at the receiver side. We consider separate detection and decoding and evaluate several options for soft-output symbol detection in terms of their pragmatic capacity, i.e., the achievable rate of the channel induced by the signal constellation used with uniform probability at the input and the detector soft output [23, 24]. This performance metric characterizes information theoretical rates achievable by separate detection and decoding, when a specific signal constellation and a specific soft-output symbol detector are employed. Therefore, it is more meaningful and practically relevant than uncoded bit error rate, as usually considered in concurrent works. We propose a message-passing (MP) detector derived from the general approach of [25], suitably adapted to the OTFS signal format. The proposed scheme is compared with the following schemes: i) another MP-based scheme recently proposed in [19]; ii) the standard linear minimum mean square error (MMSE) block equalizer, which is prohibitively complex due to its computation of a large-dimensional matrix inversion; iii) a recently proposed low-complexity approximated linear MMSE equalizer [26], based on some drastic simplifying assumptions that provide a conveniently structured channel matrix for which the large dimensional matrix inversion can be avoided. Our simulation results show that the proposed detection scheme significantly outperforms all other low-complexity schemes, and it is also able to outperform the high-complexity MMSE block equalizer when the channel is sufficiently sparse.

The main contributions of this work are summarized as follows:
1) we propose an efficient algorithm to estimate range and velocity for OTFS modulation employing practical rectangular pulses over a -path time-frequency selective radar channel. The proposed algorithm is different from the one recently proposed in [14]. While [14] focuses on a low-complexity matched filter approach, we consider ML parameter estimation. Our proposed algorithm coincides with the ML estimator for a single path channel () and it performs also very well for . Furthermore, our proposed algorithm is able to handle continuous values of delay and Doppler, contrary to [14] that assumes discretized delays.
2) We derive the CRLB of range and velocity estimation for a general -path time-frequency selective channel. For the special case of a single-path channel (), we also analyze the typical transition behavior of the ML estimator (usually referred to as waterfall), following an approach similar to [27, 28].
3) Through these analytical results and of computer simulations, we show that fully digitally modulated waveforms such asOTFS and OFDM provide as accurate range and velocity estimation as FMCW, one of typical radar waveforms. Furthermore, OTFS achieves slightly better data rates than OFDM since the latter incurs a higher overhead due to the cyclic prefix (CP).
4) We adapt the MP-based detection scheme of [25] to the OTFS format, evaluating its performance in terms of pragmatic capacity. Our scheme improves upon existing works in various aspects. First, we consider OTFS modulation with practical rectangular pulses, rather than ideal pulses satisfying the bi-orthogonal condition as considered in [18, 26, 29]. Unfortunately such ideal pulses are mathematically impossible to construct [30]. Second, our method constructs a (e.g. [31]) of girth 6, for which the exact sum-product algorithm (SPA) (e.g. [31]) can be directly applied with linear complexity in the symbol constellation size, as opposed to the MP detector of [19], based on a FG of girth 4, and requiring a Gaussian approximation of the interfering constellation symbols in order to avoid high complexity computations (more details on this point are discussed in Section V). Finally, our scheme does not make simplifying assumptions on the channel model such as discretized delay and Doppler shifts and/or exactly bi-orthogonal pulses, therefore it is much more robust to realistic system and channel conditions.

The paper is organized as follows. In Section II, we present the physical model. In Section III, we derive the OTFS input-output relation and the definition of the channel matrix. In Section IV, we provide the ML radar estimator together with the theoretical analysis of the ML estimation performance. Section V introduces the considered symbol detection algorithms. In Section VI, we define the chosen performance metrics and then we show some illustrative numerical results. Finally, Section VII concludes the paper.

## Ii Physical model

We consider a joint radar and communication system over a channel bandwidth operating at the carrier frequency . We assume that a transmitter, equipped with a mono-static full-duplex radar, wishes to convey a message to its target receiver while estimating parameters of interest related to the same receiver from the backscattered signal. Full-duplex operations can be achieved with sufficient isolation between the transmitter and the (radar) detector and possibly interference analog pre-cancellation in order to prevent the (radar) detector saturation [32]. For simplicity, in this paper we consider no self-interference, keeping in mind that some residual self-interference can be handle as additional noise [12].

We model the channel as a -taps time-frequency selective channel given by

 h(t,τ)=P−1∑p=0hpej2πνptδ(τ−τp), (1)

where denotes the number of scattering reflections or paths and the 0-th path is LoS, is a complex channel gain including the path loss (PL) of the path component, while , , , and denote the corresponding one-way Doppler shift, delay, velocity, and range associated to the -th propagation path, respectively. The round-trip Doppler shift and delay at the radar detector (co-located with the transmitter) are obtained by doubling the one-way values of the LoS path.

We assume that the transmitter periodically scans narrow angular sectors (e.g., using a phased array). If the sector beam is narrow enough, it is reasonable to consider a single target in each sector, while the propagation can go through multiple paths, one of which is the LoS, and the others may be caused by ground reflection and reflections on buildings and metal surfaces (e.g., vehicles parked along the street). An example of this scenario is depicted in Fig. 1. Therefore, we focus on the estimation of the delay (related to the range between transmitter and receiver) and of the Doppler shift (related to the relative velocity between transmitter and receiver) of the LoS path, since in the envisaged scenario the direction of arrival is automatically determined by the scanned angular sector.

## Iii OTFS Input-Output Relation

We derive here the input-output discrete complex baseband model for the OTFS signal format by generalizing previous works [18, 19], removing certain oversimplifying assumptions such as bi-orthogonal pulses and discretized delay and/or Doppler shifts.

As usually done in OTFS modulation (e.g., [18, 19]), data symbols , for and , are arranged in an two-dimensional grid referred to as the Doppler-delay domain. We consider the average power constraint given by

 1NMN−1∑k=0M−1∑l=0E[|xk,l|2]≤Pavg. (2)

In order to send the block of symbols , the transmitter first applies the inverse symplectic finite Fourier transform (ISFFT) to convert data symbols into a block of samples in the dual domain, referred to as the time-frequency domain. This is given by

 X[n,m]=N−1∑k=0M−1∑l=0xk,lej2π(nkN−mlM), (3)

for and . Then, it generates the continuous-time signal

 s(t)=N−1∑n=0M−1∑m=0X[n,m]gtx(t−nT)ej2πmΔf(t−nT), (4)

where denotes the transmit shaping pulse, is the symbol duration, and we consider a multi-carrier system where the total bandwidth is divided into subcarriers, i.e., , where denotes the subcarrier spacing. By considering strictly time limited in , the OTFS frame duration is given by (generally the duration is larger, but in practice it is well approximated by up to some guard interval containing the “tail” of ). The noiseless received signal after transmission through the time-frequency selective channel in (1) is given by

 r(t)=∫h(t,τ)s(t−τ)dτ=P−1∑p=0hps(t−τp)ej2πνpt, (5)

 Y(t,f)=∫r(t′)g∗rx(t′−t)e−j2πft′dt′. (6)

By sampling at and , the received samples in the time-frequency domain are given by

 Y[n,m]=Y(t,f)|t=nT,f=mΔf=N−1∑n′=0M−1∑m′=0X[n′,m′]Hn,m[n′,m′], (7)

where, by letting , we have

 Hn,m[n′,m′]≜P−1∑p=0h′pej2πn′Tνpe−j2πmΔfτpCgtx,grx((n−n′)T−τp,(m−m′)Δf−νp), (8)

where denotes the cross-ambiguity function between two generic pulses and , as defined in [30]. Finally, the received samples in the Doppler-delay domain are obtained by applying the symplectic finite Fourier transform (SFFT) to (7), i.e.,

 y[k,l]=1NMN−1∑n=0M−1∑m=0Y[n,m]e−j2π(nkN−mlM)=N−1∑k′=0M−1∑l′=0xk′,l′gk,k′[l,l′], (9)

where the inter-symbol interference (ISI) coefficient of the Doppler-delay pair seen by sample is given by

 gk,k′[l,l′]=P−1∑p=0h′pΨpk,k′[l,l′], (10)

where the channel matrix entries are defined in (III).

Writing the matrices of transmitted symbols and received samples as

-dimensional column vectors (stacking the columns of the corresponding matrices on top of each other), we obtain the block-wise input-output relation as

 y=(P−1∑p=0h′pΨp)Ψx+w, (12)

where is the matrix obtained from (III) while denotes the additive white Gaussian noise (AWGN) with zero mean and covariance . Notice that our input-output relation in (12) is exact (i.e., no approximation was made in its derivation) and holds for any pair of transmit/receive pulses.

In order to proceed further in a tractable manner, as done in [19], we approximate the integral of the cross-ambiguity function with a discrete sum. In particular, we write

 Cgrx,gtx(τ,ν) =∫T0gtx(t)g∗rx(t−τ)e−j2πνtdt ≈TMM−1∑i=0gtx(iTM)g∗rx(iTM−τ)e−j2πνiTM, (13)

where, for analytical convenience, the number of equally spaced discretization nodes is equal to the number of subcarriers . This approximation is accurate provided that large enough, which is the typical case. Now, by letting and be rectangular pulses of length , by limiting the channel delay to , it readily follows that the cross-ambiguity function yields non-zero samples only for and .

Hence, the entries of the cross-talk matrix can be thus written as in (III), in which we exploited the expression of the approximated cross-ambiguity function for rectangular pulses given by

 Cgrx,gtx(τ,ν)≈1MM−1−lτ∑i=0exp(j2πνiTM), (15)

with , and where we defined the index sets

 {LICI≜[0,M−1−lτp]LISI≜[M−lτp,M−1]. (16)

Moreover, in order to have the received signal bandwidth still approximately equal to , we must assume that the bandwidth expansion due to Doppler is negligible with respect to . In our case, for the sake of simplicity, we make the usual assumption that .111Note that this approximation can be justified in a number of scenarios. For example, consider a scenario inspired by IEEE 802.11p with GHz and the subcarrier spacing KHz. This yields km/h, which is reasonable even for a relative speed of 400 km/h. The same holds for IEEE 802.11ad with = 60 GHz and MHz [33].

### Iv-a Maximum Likelihood Estimator

Focusing on the channel model defined in (1) with a single target and multiple backscattered paths, we wish to find the ML estimator for the set of unknown parameters , where the bar indicates the true value of the unknown parameters.

By using expression (III) for the channel matrix coefficients, which contains the dependency on the parameters , the log-likelihood function to be minimized is given by

 l(y|θ,x)=∥∥ ∥∥y−P−1∑p=0h′pΨpx∥∥ ∥∥2, (17)

where symbols in are known at the radar detector, since it is co-located with the transmitter. The ML estimator is given by

 ^θ=argminθ∈CP×RP×RPl(y|θ,x). (18)

A brute-force search for the maximum in a dimensional continuous domain is infeasible in general. Therefore, in the following we propose a viable method to approximate the ML solution with low complexity.

The log-likelihood function in (17) is quadratic in the complex amplitudes for given . Hence, the minimization of (17) with respect to for fixed is readily given as the solution of the linear system of equations

 P−1∑q=0h′qxHΨHpΨqx =xHΨHpy,p=0,…,P−1. (19)

Expanding (17) and using the equality (19), after some long but relatively simple algebra (not given explicitly for the sake of brevity), we find that the minimization with respect to reduces to maximizing the function

 l2(y| θ,x)=P−1∑p=0⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩|xHΨHpy|2xHΨHpΨpxSp(τp,νp)−(∑q≠ph′qxHΨHpΨqx)yHΨpxxHΨHpΨpxIp({h′q}q≠p,τ,ν)⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭, (20)

where and denote the useful signal and the interference for path , respectively. Clearly, since the channel coefficients are not known, it is impossible to directly maximize w.r.t. . Furthermore, even for known coefficients , the function is not separable in the pairs of parameters for different values of because of the dependency of the interference terms on all for . Nevertheless, this dependency appears through the “cross-term” coefficients of the type , which tend to be weak for typically sparse multipath channels. Therefore, we resort to an iterative block-wise optimization that alternates the optimization of each pair by keeping fixed the other parameter pairs and the channel complex coefficients and, after a round of updates for the delay and Doppler parameters, it re-evaluates the estimates of the channel coefficients by solving (19). The algorithm steps are given in the following, where the iteration index is denoted in brackets by :

1. Initialization: let and for all .

2. For repeat:

• Delay and Doppler shift update: for each , find the estimates by solving the two-dimensional maximization

 (^τp[nit],^νp[nit]) = argmax(τp,νp){Sp(τp,νp)
• Complex channel coefficients update: solve the linear system (19) for the channel matrices calculated with parameters , and let the solution be denoted by .

The iteration stops when the values of do not show a significant change with respect to the previous iteration, or if a maximum number of iterations is reached. In practice, we find the maximizer in (LABEL:diocane) by searching on a finely discretized grid on the delay and Doppler domain. Furthermore, in all our simulations we noticed that the algorithm converges in just a few iterations (2 to 5, at most).

### Iv-B Performance Bounds

#### Iv-B1 Cramér-Rao lower bound (Crlb)

The derivation of the CRLB is based on the channel matrix expression in (III). From the channel model (1), by letting

 sp[k,l]=∣∣h′p∣∣ej∠h′pN−1∑k′=0M−1∑l′=0Ψpk,k′[l,l′]xk′,l, (22)

and considering separately the amplitude and the phase of the complex channel coefficients , the element of the Fisher information matrix is given by

 [I(θ,x)]i,j =2Pavgσ2wRe{∑n,m[∂sp[n,m]∂θi]∗[∂sp[n,m]∂θj]}, (23)

in which is the set of channel unknown parameters.

The partial derivatives w.r.t. the magnitude and phase of are straightforward and are omitted for the sake of space limitation. The derivatives w.r.t. and w.r.t. are more cumbersome and, after some algebra, can be obtained in expressions (IV-B1) and (IV-B1), respectively.

The desired CRLB follows by filling the Fisher information matrix in (23) with the derivatives computed above, and obtaining the diagonal elements of the inverse Fisher information matrix. In particular, we are interested in the CRLB for the parameters and , related to the target range and velocity.

#### Iv-B2 Waterfall Analysis for P=1

In the case of our iterative scheme reduces to the exact joint ML estimation of (in other words, the iterative scheme converges after the first iteration to the exact maximum of the log-likelihood function, up to the discretization error in the 2-dimensional search domain ). It is well-known that ML estimators typically exhibit a threshold effect, i.e., a rapid deterioration of the estimation MSE when the signal-to-noise ratio (SNR) is below some threshold (that generally depends on the problem and on the sample size). In contrast, for SNR larger than such threshold the ML estimator yields MSE typically very close to the CRLB

. This waterfall behavior is caused by “outliers” in the search of the maximum in

ML estimator: namely, when the observation is too noisy, the maxima of the likelihood functions tend to be randomly placed anywhere on the search grid. In the following we provide an analysis of the waterfall transition for the case , i.e., the region around the threshold SNR value where the rapid deterioration occurs. Our simulations show that the waterfall prediction provided by our analysis for is also very accurate for the multipath case . This corroborates the evidence that our proposed approximated ML algorithm is effectively very good, and performs very close to the true ML.222Obviously, the CRLB for yields also a lower bound for the case , in the case where we simply add more multipath components on top of the LoS component without changing the statistics of the LoS component, since the presence of more unknown parameters cannot help the estimation of the target parameters .

Following the reasoning of [27, 28], we treat as “outlier” the event that a maximum of the log-likelihood function is randomly placed on the grid , rather than in the cluster gird points around the true value . Let

be the unknown parameter to be estimated. By the law of total probability over the discretized grid

(where we calculate the ML estimator), we can approximate the estimation MSE as

 MSE=E[(^α−¯α)2] =∑i∈ΓPr(ϵ(i))(^αi−¯α)2≤∑i∈ΓPr(~ϵ(i))(^αi−¯α)2. (26)

where is the true value of the parameter, is the estimated parameter, and denotes the probability of error of choosing rather than . While evaluating may be extremely difficult, we obtain an upper bound by considering pairwise error probabilities, i.e., replacing with the probability that the detector chooses rather than when these are the only two alternatives. Since the pairwise error event contains the true error event , it follows that the inequality of (26) provides an upper bound.

Proceeding further, we define the pairwise error probability as

 Pr(~ϵ(i))≜Pr{l(y|θi,x)

where and are the value obtained from the evaluation of the likelihood function at grid points with parameters and using the true parameters , respectively. At this point, the problem is reduced to the computation of the pairwise error probabilities , for , which can be derived as follows. We notice that

 Pr(~ϵ(i))≈Pr{|xHΨHiy|2>|xH¯ΨHy|2}, (28)

where is the channel matrix calculated with the true parameters and we exploit , , for

. Defining the jointly conditionally Gaussian random variables

 zi ≜ xHΨHiy=xHΨHi¯Ψx+xHΨHiw (29) ¯z ≜ xH¯ΨHy≈xHx+xH¯ΨHw, (30)

with first and second order moments

 {E[¯z]=∥x∥2,Var[¯z]=σ2w∥x∥2E[zi]=xHΨHi¯Ψx,Var[zi]=σ2w∥x∥2, (31)

and

 Cov[¯z,zi]=σ2wxH¯ΨHΨix, (32)

a good approximation for is given by [27]

 Pr(~ϵ(i))≈12exp{−∥x∥42σ2wNM}I0(∣∣xHΨHi¯Ψx∣∣⋅∥x∥22σ2wNM), (33)

where is the modified Bessel function of the first kind of order . The MSE estimation can be thus computed by substituting (33) into (26), using the above definitions.

The resulting (approximated) upper bound tends to be loose at low SNR, but becomes more accurate moving towards the ML waterfall region. Furthermore, the MSE strictly depends on the grid resolution and may or may not reach the CRLB for increasing SNR depending on the systematic error incurred by the grid discretization. As a matter of fact, in our numerical results we took care of using a search grid fine enough such that the discretization systematic error is not visible in the explored range of SNR.

In order to cope with the fact that the (approximated) MSE upper bound obtained in (26) becomes very loose for low SNR, we use also the trivial upper bound

 MSE≤1|Γ|∑i∈Γ(^αi−¯α)2, (34)

that corresponds to choosing at random grid points, by disregarding completely the received signal.

Then, putting together (26) (with the pairwise error probability approximation in (33) and the above “random estimation” bound in (34), we finally obtain the approximated MSE upper bound

 (35)

Our simulations show that this approximated upper bound is very accurate and it is able to predict very well the waterfall behavior of the ML estimator (as well as the proposed approximated ML estimator in the case ).

## V OTFS Soft-Output Data Detection

We will now focus on the OTFS data detection at the receiver side. By considering the communication channel model in (1), in line with most of the current literature on OTFS detection (e.g., [19, 18]), we assume perfect channel state information (CSI) at the receiver. An efficient pilot-aided CSI acquisition is a problem of independent interest that we postpone to future work (see also [21, 20, 34]).

In this paper, we consider separate detection and decoding, where the receiver consists of the concatenation of a soft-output symbol detector, producing soft-estimates of the coded symbols , and a decoder that takes such estimates as the output of a virtual channel that incorporates also the detector. Motivated by practical complexity considerations, no “turbo equalization” involving a feedback loop from the (soft-output) decoder to the detector and iterations between detector and decoder is considered. It follows that the relevant performance measure is the already mentioned pragmatic capacity, i.e., the mutual information between the input constellation symbols, used with uniform probability, and the corresponding detector soft-output [23, 24].

We propose an efficient low-complexity MP-based soft-output detector, obtained by constructing a FG

for the joint posterior probability of

given in (12) and applying the standard SPA computation rules to compute its marginals. The FG is constructed according to the general approach of [25] (applicable to any linear-Gaussian model such as (12)), for which the graph girth is guaranteed to be at least 6 and the degree of the function nodes is at most 2. This allows the application of the exact SPA computation at the nodes and a high degree of parallelization, such that the resulting MP-based detector is very computationally efficient. Furthermore the absence of cycles of length less than 6 yields good convergence properties of the SPA iterations.

We compare the proposed scheme with other three soft-output detectors proposed in the literature, namely: i) a different MP-based detector recently proposed in [19], based on an alternative way to construct the FG (see Section V-B); ii) a standard linear MMSE block equalizer, which offers very good performance at the impractical cost of a very large-dimensional matrix inversion; iii) a low-complexity MMSE block equalizer recently proposed in [26] that uses drastic simplifying assumptions (in particular, the bi-orthogonality of the OTFS and pulses and the assumption that the delay and Doppler shifts as integer multiples of the receiver sampling grid) such that the resulting nominal channel matrices are block-circulant with circulant blocks and the matrix inversion in linear MMSE estimation can be efficiently implemented. As we shall see, since these simplifying assumptions are not verified in practice (shifts are not on a quantized grid, and bi-orthogonal pulses with unit time-frequency product cannot exist [30]), the performance of this low-complexity linear detector is quite poor, when applied to a realistic channel and pulse scenario.

### V-a Proposed MP-based detector (“Matrix G algorithm” — MPG)

From (12), when is known (perfect CSI) and is complex Gaussian i.i.d., the conditional probability density function (pdf) of the received samples given the modulation symbols is given by

 (36)

where the proportionality symbol indicates an irrelevant constant factor independent of symbols . We follow the FG construction approach of [25], and expand the -norm inside the exponential as

 ∥y−Ψx∥2=yHy−2Re{xHΨHy}+xHΨHΨx. (37)

Defining and , the conditional pdf can be written as

 p(y|x)∝exp(2Re{xHz}−xHGz2σ2w). (38)

Note that the sequence is a sufficient statistic for symbols detection. By expressing the matrix operations explicitly in terms of the components, we define the functions

 Fi(xi)=exp[1σ2wRe{zix∗i−Gi,i2|xi|2}], (39)
 Ii,j(xi,xj)=exp[−1σ2wRe{Gi,jxjx∗i}], (40)

and we use Bayes rule in order to express the a-posteriori probability of given in the factored form, i.e.,

 (41)

where we used the fact that the modulation symbols take on values in some signal constellation and are treated by the detector as i.i.d. with given (typically uniform) a priori probability mass function .

In the proposed approach, the FG corresponds to the factorization in (41) (see the example shown in Fig. 2). At this point, the resulting MP-based soft-output detector follows immediately by applying the standard SPA computation rules. The detailed derivation of the message computation at the function and variable nodes of the FG is given in [25] and applies directly to our setting. Here, for the sake of completeness, we just summarize the resulting algorithm. Define as the product of all messages incoming to the variable node , namely

 Vi(xi)=Pi(xi)Fi(xi)∏j≠iνi,j(xj), (42)

which is proportional to the (estimated) a-posteriori probability and thus provides the soft-output of the this detector. Then, the application of the SPA leads to the following rules for message exchange and update:

1. Computation at the variable nodes: each node sends to each adjacent function node the message

 μi,j(xi)=Vi(xi)/νi,j(xi). (43)
2. Computation at the function nodes: each function node sends to each adjacent variable node the message

 νi,j(xi)=∑xj∈CIi,j(xi,xj)μi,j(xj). (44)

Notice that all the variable nodes and the function nodes can be activated in alternative rounds and, in each round, all the nodes of the same type can be activated in parallel, e.g., adopting the same flooding schedule used for low-density parity check (LDPC) decoding [35]. Moreover, messages can be implemented in the logarithmic domain [31].

With the help of Fig. 2, we can illustrate some important features of the proposed approach.

• FG of girth 6: the FG constructed as above is guaranteed to have girth (minimum length of cycles) equal to 6 (highlighted in bold in Fig. 2). It is well-known that the SPA yields exact posterior marginalization for cycle-free FGs, and the rationale behind the use of the SPA paradigm on loopy graphs is that if the FG has large girth, the local neighborhood of each node is “tree-like” [31]. In particular, cycles of length 4 should be avoided. Hence, the proposed construction following the general method of [25] yields indeed an FG better suited to the application of the iterative SPA.

• Computational complexity: notice that the computation in (40) involves the summation with respect to only a single discrete variable over the constellation . Therefore, this computation has always linear complexity in the constellation size, irrespectively on the sparsity of the channel matrix . This allows the use of the exact SPA at fixed complexity per node, unlike the approach in [19] (see comments in Section V-B).

• High degree of parallelization: the number of nodes of type depends on the number of non-zero elements in the rows of the upper triangular part of the matrix (the existence of edges is evidenced by the shadowed elliptic area in Fig. 2). Nevertheless, as said before, these degree-2 nodes can be all activated in parallel. Hence, for a sufficiently large degree of parallelization, the computational time complexity is independent of the sparsity of the multipath channel. Notice that in modern LDPC decoding is not unlikely to find implementations with degree of parallelization of the order of 1000, which is much larger than what needed in our detector. Hence, we claim that the proposed detector is very attractive from a practical implementation viewpoint.

### V-B MP-based algorithm of [19] (“Matrix Ψ algorithm” — MPΨ)

The MP-based algorithm of [19] builds its FG from the more “direct” factorization of the a-posteriori probability

 (45)

Defining the function nodes

 Qi(x)≜p(yi|x)=exp(∥yi−Ψix∥22σ2w), (46)

where is the -th row of , the resulting FG is shown in Fig. 3. Since the elements of the same -th column is multiplied by the same symbol , this FG has necessarily cycles of length 4 (highlighted in bold in Fig. 3). The edges evidenced by the shadowed elliptic area in the figure correspond to the non-zero elements of the matrix . In particular, the degree of the function node is equal to the number of non-zero elements in the the -th row . The exact SPA computation at such function nodes requires summing over discrete variables taking values in . Therefore, it has complexity that may be prohibitively large for large constellations and, above all, it depends on the sparsity of the channel. Therefore, the exact application of the SPA computation rules to the FG obtained directly from the matrix is highly impractical. For this reason, the authors of [19] propose to use a Gaussian approximation of the interfering symbols in the computation at the nodes , which effectively boils down to a soft interference cancellation approach, as already widely used in turbo equalization and pioneered in the context of multiuser detection in [36, 37]. For the sake of space limitation, we omit the details of the resulting MP algorithm, which can be found in [19].

It should also be mentioned that, for the sake of simplicity and in order to increase the sparsity of the FG in Fig. 3, the detector proposed in [19] constructs the nominal matrix by rounding the delay shifts to integers on receiver sampling grid. Under this condition, the channel matrix is very sparse since many coefficients corresponding to sampling at non-integer delay shifts are identically to zero, and the number of connections for each node is reduced, while preserving the 4-cycle problem said before since this is unavoidable with this approach. Nevertheless, since the assumption is generally not satisfied by real-world channels, such approximation of the channel matrix results in neglecting a significant component of the ISI. We shall verify that when such integer delay shift rounding is applied to the construction of the nominal matrix used by the detector, but the actual delay shifts have a fractional component (as it is always the case in practice), this mismatch yields a significant performance degradation. This shows that neglecting the fractional part of delays and Doppler shifts, as routinely done in the literature of OTFS, may be indeed quite misleading.

### V-C Linear block-wise MMSE equalization

As a further term of comparison we consider also the standard linear MMSE block equalizer, applied to the channel model (12). In this case, the soft-output is simply the linear MMSE estimate of symbols from the observation , given by

 ^xLMMSE=ΨH(ΨΨH+σ2wI)−1y. (47)

The complexity of this approach is proportional to , so it becomes quickly unfeasible for typical values of and (e.g., in our simulations we have considered and ). A low-complexity (mismatched) linear minimum mean square error (LMMSE) approach was recently proposed in [26]. This relies on the cyclic properties of channel matrix under perfect bi-orthogonality of the modulation pulses and and integer-grid valued delays and Doppler shifts. As a result, in the presence of practical rectangular pulses non-integer delays and Doppler shifts, as realistically considered in our work, the performance of this approach visibly degrades.

## Vi Simulation Results

The radar (backscattered) and forward communication channels are both defined by (1) with different values of PL, Doppler shift, and delay. We define the radar and communication SNRs as [38, Chapter 2]

respectively, where is the wavelength, is the radar cross-section in , is the antenna gain, and is the distance between transmitter and receiver.

In the case of multipath, we fix to be the SNR of the LoS component, and we add multipath components with progressively lower SNRs, such that the sum SNR of the channel increases with the number of paths . This corresponds to the physically meaningful case that a richer propagation environment conveys more signal power. Table I summarizes the relevant simulation parameters inspired by the automotive communication standard IEEE 802.11p [6], where and denote the target range and velocity.

### Vi-a Comparison with OFDM and FMCW

We briefly review OFDM and the widely used radar waveform known as FMCW [38, Chapter 4.6]. For both OFDM and FMCW we consider a symbol length of including a guard interval denoted by , longer than the maximum path delay . In OFDM, the CP length is given by , with . We send modulation symbols in the time-frequency domain, satisfying the average power constraint (2). The OFDM receiver samples at rate and removes the CP before detection in order to eliminate the inter-symbol interference. For the case of a single path channel (), the ML estimator and the related CRLB can be found in our related work [22] and are omitted here for the sake of space limitation.

In FMCW, the radar transmitter sends a sequence of identical “chirp” pulses of duration , each followed by a guard interval of to avoid inter-pulse interference. By letting denotes the phase at time , the transmit signal is given by

 s(t)=N−1∑i=0ej2πϕ(t−iT0)rect(t−iT0T), (49)

where we consider pulses to make a fair comparison with OFDM. Plugging (49) into (5) , we obtain the received signal by neglecting the noise. After some algebra, it is easy to show that the product of the received and the transmit signals gives

 y(t) =r(t)s∗(t)=P−1∑p=0hpej2πνpte−j2πfcτpN−1∑i=0e−j2πfb,p(t−iT0−τp/2), (50)

where we let denotes the so-called beating frequency of path . The receiver samples every for each pulse, i.e. for where denotes the pulse index and denotes the sample index. By letting the number of samples par pulse, the sampled received signal can be rewritten as

 y[i,l] Δ=y(iT0+lfs)=∑phpej2π(fb,p+νp)lfsej2πνpiT0, (51)

for and , where absorbs a constant phase term independent of the indices .

The estimation of the unknown parameters is thus obtained by selecting the peak of the range-Doppler map found by applying a two-dimensional discrete Fourier transform (DFT) to the noisy samples in (51) as proposed in [12, 4]. The range and velocity of the target are obtained by the estimates of .

### Vi-B Joint Radar and Communication Performance

The first two subfigures of Fig. 4 show the velocity and range estimation root MSE (RMSE) versus for a pure LoS channel () and for OTFS, OFDM, and FMCW. We notice that both digital modulation formats provide as accurate radar performance as FMCW, while transmitting at their full information rate. It is remarkable to note that the combination of the waterfall analysis and the CRLB introduced in Section IV-B are able to accurately predict the actual ML estimation performance.

In addition, the third subfigure of Fig. 4 shows the achievable rate with Gaussian independent and identically distributed inputs symbols