# Monostatic Sensing with OFDM under Phase Noise: From Mitigation to
Exploitation

Musa Furkan Keskin, Member, IEEE, Henk Wymeersch, Senior Member, IEEE, and
Visa Koivunen, Fellow, IEEE Musa Furkan Keskin and Henk Wymeersch are with the
Department of Electrical Engineering, Chalmers University of Technology, SE
41296 Gothenburg, Sweden (e-mail: furkan@chalmers.se). Visa Koivunen is with
the Department of Signal Processing and Acoustics, Aalto University, FI 00076
Aalto, Finland. This work is supported, in part, by MSCA-IF grant 888913
(OTFS-RADCOM) and the European Commission through the H2020 project Hexa-X
(Grant Agreement no. 101015956).

###### Abstract

We consider the problem of monostatic radar sensing with orthogonal frequency-
division multiplexing (OFDM) joint radar-communications (JRC) systems in the
presence of phase noise (PN) caused by oscillator imperfections. We begin by
providing a rigorous statistical characterization of PN in the radar receiver
over multiple OFDM symbols for free-running oscillators (FROs) and phase-
locked loops (PLLs). Based on the delay-dependent PN covariance matrix, we
derive the hybrid maximum-likelihood (ML)/maximum a-posteriori (MAP) estimator
of the deterministic delay-Doppler parameters and the random PN, resulting in
a challenging high-dimensional nonlinear optimization problem. To circumvent
the nonlinearity of PN, we then develop an iterated small angle approximation
(ISAA) algorithm that progressively refines delay-Doppler-PN estimates via
closed-form updates of PN as a function of delay-Doppler at each iteration.
Moreover, unlike existing approaches where PN is considered to be purely an
impairment that has to be mitigated, we propose to exploit PN for resolving
range ambiguity by capitalizing on its delay-dependent statistics (i.e., the
range correlation effect), through the formulation of a parametric Toeplitz-
block Toeplitz covariance matrix reconstruction problem. Simulation results
indicate quick convergence of ISAA to the hybrid Cramér-Rao bound (CRB), as
well as its remarkable performance gains over state-of-the-art benchmarks, for
both FROs and PLLs under various operating conditions, while showing that the
detrimental effect of PN can be turned into an advantage for sensing.

Index Terms– OFDM, joint radar-communications, phase noise, exploitation,
iterated small angle approximation.

##  I Introduction

Envisioned as one of the key enabling technologies in 6G wireless networks,
the concept of joint radar-communications (JRC), also known as integrated
sensing and communications (ISAC), has drawn significant attention in recent
years [1, 2, 3, 4, 5, 6]. Towards practical JRC implementation in real-world
scenarios, two principal design paradigms have been pursued in the literature:
radar-communications coexistence (RCC) [7] and dual-functional radar-
communications (DFRC) [8, 9]. RCC deals with joint optimization of spectrally
coexisting radar and communication systems deployed on separate platforms to
mitigate mutual interference and enable efficient spectrum utilization [10],
while DFRC refers to integration of sensing and communication functionalities
into a single hardware that employs a joint waveform to perform both tasks
simultaneously for mutual benefit [4]. Sparking innovative use cases (e.g.,
sensing-assisted communications [6]) and bringing hardware/cost efficiency by
piggybacking on the existing wireless infrastructure, DFRC holds great
potential in emerging 6G networks, where sensing will be an integral component
[6, 11]. As a promising candidate for DFRC transmission, the orthogonal
frequency-division multiplexing (OFDM) waveform has been extensively studied
thanks to its wide prevalence in mobile network standards and its satisfactory
performance in radar operation [12, 13, 14, 15].

As next-generation wireless systems are expected to operate at high
frequencies, i.e., millimeter wave (mmWave) bands [11], hardware impairments
(HWIs), such as phase noise (PN) [16], power amplifier nonlinearity (PAN) and
mutual coupling (MC), can become a major bottleneck for OFDM DFRC system
performance in both radar and communications [17, 18]. In particular, the
severity of PN, which is caused by oscillator imperfections111Being a time-
varying impairment, PN constitutes a much more serious issue for DFRC systems
than static impairments caused by oscillator non-idealities, such as carrier
frequency offset (CFO) and I/Q imbalance [19]. For instance, CFO has no effect
on monostatic sensing systems as the same oscillator is employed for
transmission and reception [20]., increases with operating frequency [19, 21].
Due its rapidly time-varying nature, PN requires dynamic compensation
accounting for subcarrier-level (frequency-domain) or sample-level (time-
domain) mitigation processing, thereby posing a significant challenge in
channel estimation and data detection for OFDM communications [22, 19, 23,
24]. To tackle the PN compensation problem for OFDM, various frequency-domain
[25, 21, 26, 27, 28] and time-domain [29, 30, 16] estimation approaches have
been proposed.

Despite the vast literature on PN estimation in OFDM communications, very
little effort has been devoted to studying the impact of PN on the performance
of OFDM radar (e.g., [31, 32]), let alone to developing algorithms to estimate
and compensate for PN in radar sensing. In [31] and [32], the effect of PN on
the range-velocity profile of an OFDM radar is investigated, which shows that
PN leads to an increase in noise floor and produces a ridge along the velocity
axis. To support new use cases towards 6G networks, emerging mmWave sensing
applications [33, 34, 35, 36] impose stringent requirements on range and
velocity accuracy [37], which necessitates the consideration of PN. From the
perspective of radar receive processing, the existing OFDM radar algorithms
(e.g., [12, 38, 39, 40, 14, 41]) assume ideal oscillators and thus cannot
provide satisfactory performance under the effect of PN, especially for large
PN variances. In a nutshell, no systematic study has been performed to tackle
the problem of radar sensing in OFDM DFRC systems with oscillator PN, and,
accordingly, to derive an algorithm to jointly estimate delay, Doppler and PN.

As a key distinction between radar and communications, the so-called range
correlation effect [42, 43, 44, 45] constitutes an essential peculiarity of PN
in monostatic radar sensing compared to PN in a communications setup. Due to
independent PN processes at the transmitter and the receiver, the PN
statistics in communications systems do not depend on unknown channel
parameters (e.g., [46, 47, 19, 16]). Conversely, in a shared-oscillator radar
transceiver, downconversion of the reflected signal from a target results in a
differential (self-referenced/self-correlated) PN process [48, 49]
corresponding to the difference between the original PN process and its
version shifted in time by the round-trip delay of the target. This
correlation effect renders the statistics of PN in the radar receiver range-
dependent, leading to higher PN variance for farther targets [43], which
brings both challenges and opportunities specific to radar sensing. The main
challenge pertains to delay and PN being coupled, making it difficult to
disentangle the corresponding estimation tasks (as is commonly performed in
joint channel/PN estimation in OFDM communications, e.g., [29, 50]). On the
other hand, the main opportunity arises from the possibility to exploit the
range-dependent PN statistics for enhancing range estimation performance.
While PN exploitation may offer promising performance gains, to the best of
authors’ knowledge, this topic remains surprisingly unexplored, both in the
standard radar literature (i.e., frequency-modulated continuous wave (FMCW),
multiple-input multiple-output (MIMO) or pulsed radars) and in OFDM DFRC
research.

In light of the existing literature, several fundamental questions arise
concerning the sensing functionality of OFDM DFRC systems in the presence of
PN:

  * •

What are the statistical properties of PN in the OFDM radar receiver for
different types of oscillators, namely, free-running oscillators (FROs) and
phase-locked loops (PLLs) [48, 25]?

  * •

How can we develop powerful algorithms for joint estimation of delay, Doppler
and PN, to mitigate the impact of PN on the sensing performance? How do FROs
and PLLs affect the performance of delay and Doppler estimation?

  * •

Considering the range correlation effect, is it possible to exploit PN to
improve ranging performance beyond that achievable via PN-free, ideal
oscillators?

In an attempt to answer these questions, this paper studies the problem of
radar delay-Doppler estimation in OFDM DFRC systems in the face of oscillator
PN. We begin by deriving the statistical characteristics of PN in the OFDM
radar observations for both FROs and PLLs. Then, we propose a novel algorithm
for joint estimation of delay, Doppler and PN, based on iterated small angle
approximation of PN in the cost function of the hybrid maximum-likelihood
(ML)/MAP estimator, which enables fast convergence to the corresponding
theoretical bounds. Furthermore, we develop a PN exploitation approach that
can effectively utilize the delay-dependent PN statistics to resolve range
ambiguity. The main contributions of this paper can be summarized as follows:

  * •

Problem Formulation for OFDM Radar Sensing under PN: For the first time in the
literature, we investigate the problem of monostatic radar sensing in OFDM
DFRC systems under the effect of oscillator PN. To provide a rigorous problem
formulation, we derive an explicit statistical characterization of PN in the
OFDM radar receiver. We consider two commonly used oscillator models, namely,
FRO and PLL. The derivations reveal the block-diagonal structure of the PN
covariance matrix for the former and the more general Toeplitz-block Toeplitz
structure for the latter.

  * •

Hybrid ML/MAP Estimator via Iterated Small Angle Approximation Algorithm: We
derive the hybrid ML/MAP estimator of the deterministic delay-Doppler
parameters and the random PN over an OFDM frame with multiple symbols. The
covariance matrix of PN in the backscattered signal depends on the unknown
delay. To deal with the highly nonlinear nature of the resulting cost
function, we propose a novel iterated small angle approximation (ISAA)
approach that invokes the small PN approximation around the current PN
estimate at each iteration, progressively refining delay-Doppler-PN estimates
and minimizing the impact of residual PN. The proposed approach enables
closed-form update of PN as a function of delay-Doppler and provides
significant improvements in PN tracking accuracy through alternating
iterations.

  * •

PN Exploitation to Resolve Range Ambiguity: Relying on the key insight that PN
conveys valuable information on delay through its delay-dependent statistics,
we develop an algorithm for resolving range ambiguity, that exploits the
statistics of the PN estimates at the output of the proposed ISAA method. The
PN exploitation approach formulates the range estimation as a parametric
covariance matrix reconstruction problem by leveraging the Toeplitz-block
Toeplitz structure and is capable of yielding unambiguous range estimates
through the fact that PN covariance imposes no ambiguity in range (as opposed
to the fundamental upper limit dictated by OFDM subcarrier spacing [12, 39,
51]).

  * •

Simulation Analysis: Extensive simulations conducted under a wide variety of
operating conditions indicate that the proposed ISAA algorithm converges
quickly to the corresponding hybrid Cramér-Rao bounds (CRBs) [52, 53] on
delay-Doppler-PN estimation in few iterations and considerably outperforms the
benchmark FFT method [12, 39, 40]. The accuracy gains are more pronounced for
higher SNRs, larger $3\,\rm{dB}$ oscillator bandwidths, smaller loop
bandwidths (for PLLs) and farther targets. In addition, PLLs are found to be
more beneficial for Doppler estimation than FROs, due to the presence of slow-
time PN correlation for PLLs. Moreover, above a certain SNR level, the PN
exploitation algorithm is shown to correctly identify the true range of range-
ambiguous targets and achieve much higher ranging accuracy than the FFT method
fed with PN-free observations, thereby turning PN into an advantage for
sensing.222Notations:
$\bm{\Pi}_{\mathbf{X}}=\mathbf{X}(\mathbf{X}^{H}\mathbf{X})^{-1}\mathbf{X}^{H}$
represents the orthogonal projector onto the column space of $\mathbf{X}$ and
$\bm{\Pi}^{\perp}_{\mathbf{X}}={\bm{\mathrm{I}}}-\bm{\Pi}_{\mathbf{X}}$.
$\odot$ and $\otimes$ denote the Hadamard and Kronecker product, respectively.
${\rm{diag}}\left(\mathbf{x}\right)$ outputs a diagonal matrix with the
elements of a vector $\mathbf{x}$ on the diagonals,
${\rm{diag}}\left(\mathbf{X}\right)$ represents a diagonal matrix with the
diagonal elements of a square matrix $\mathbf{X}$ on the diagonals,
${\rm{vec}}\left(\cdot\right)$ denotes matrix vectorization operator, and
${\rm{reshape}}_{N,M}\left(\cdot\right)$ reshapes a vector into an $N\times M$
matrix.

##  II System Model and Problem Formulation

Consider an OFDM JRC system consisting of a DFRC transceiver and a
communications receiver (RX), as shown in Fig. 1. Equipped with a radar-
communications transmitter (TX) (i.e., a conventional OFDM TX) and a radar RX
on a single hardware platform, the DFRC transceiver sends data symbols to the
communications RX and simultaneously performs monostatic radar sensing using
the backscattered signals to accomplish various radar tasks (e.g., target
detection, estimation, tracking and classification) [12, 8, 2]. To enable
full-duplex operation without self-interference to the radar RX, we assume
sufficient isolation and decoupling of TX/RX antennas at the DFRC transceiver
[12, 54, 39, 55, 56, 6]. At the communications RX, conventional OFDM receive
operations (e.g., channel estimation, frequency synchronization, data
detection [24]) are performed ordinarily without any constraints from the
sensing functionality. Moreover, the oscillator of the DFRC transceiver, which
is shared between the TX and radar RX on the co-designed joint hardware
platform, is assumed to be non-ideal and impaired by PN due to imperfections
[57, 46, 25, 19, 21]. In this section, we derive OFDM transmit and radar
receive signal models in the presence of PN, and formulate the resulting OFDM
radar sensing problem. We note that the paper will focus on radar sensing
under PN while the communications RX is assumed to compensate for PN via well-
established approaches, e.g., [30, 22, 58, 16].

Figure 1: OFDM JRC system comprising a DFRC transceiver that sends data
symbols to a communications RX and simultaneously performs monostatic sensing
using the backscattered signals under the impact of PN. Downconversion of the
radar return signal using the shared oscillator results in a self-correlated
PN process, leading to delay-dependent PN statistics in radar sensing, while
PN statistics in the communications RX, which employs an independent
oscillator, have no relation to channel parameters.

###  II-A Transmit Signal Model

We consider an OFDM communication frame with $M$ symbols and $N$ subcarriers.
The total duration of a symbol is given by $T_{\rm{sym}}=T_{\rm{cp}}+T$, where
$T_{\rm{cp}}$ and $T$ denote, respectively, the cyclic prefix (CP) and the
elementary symbol durations [12]. In addition, $\Delta f=1/T$ is the
subcarrier spacing, leading to a total bandwidth of $N\Delta f=B$. The complex
baseband OFDM transmit signal can be expressed as [13]

$\displaystyle s(t)=\sum_{m=0}^{M-1}s_{m}(t)\leavevmode\nobreak\ ,$ (1)

where

$s_{m}(t)=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x_{n,m}\,e^{j2\pi n\Delta
ft}{\rm{rect}}\left(\frac{t-mT_{\rm{sym}}}{T_{\rm{sym}}}\right)$ (2)

is the OFDM signal for the ${m^{\rm{th}}}$ symbol, $x_{n,m}$ denotes the
complex data symbol on the ${n^{\rm{th}}}$ subcarrier for the ${m^{\rm{th}}}$
symbol, and ${\rm{rect}}\left(t\right)$ is a rectangular pulse that takes the
value $1$ for $t\in\left[0,1\right]$ and $0$ otherwise. In the presence of PN
in the oscillator, the upconverted transmit signal over the block of $M$
symbols for $t\in\left[0,MT_{\rm{sym}}\right]$ can be written as [49]

$\widetilde{s}(t)=\Re\left\\{s(t)e^{j\left[2\pi
f_{c}t+\phi(t)\right]}\right\\}\leavevmode\nobreak\ ,$ (3)

where $f_{c}$ is the carrier frequency and $\phi(t)$ denotes the PN process in
the oscillator.

###  II-B Receive Signal Model

In radar sensing, we assume the existence of a point target in the far-field,
with round-trip delay $\tau=2R/c$, normalized Doppler shift $\nu=2v/c$ and
complex channel gain $\alpha$, which includes path loss and radar cross
section effects. Here, $R$, $v$ and $c$ denote the distance, radial velocity
and speed of propagation, respectively. Given the transmit signal model in
(3), the passband backscattered signal at the radar RX can be given as

$\displaystyle\widetilde{y}(t)=\Re\left\\{\alpha\,s(t-\tau(t))e^{j\left[2\pi
f_{c}(t-\tau(t))+\phi(t-\tau(t))\right]}\right\\}\leavevmode\nobreak\ ,$ (4)

where $\tau(t)=\tau-\nu t$ is the time-varying delay due to Doppler shift.
After downconverting the passband signal in (4) through the noisy oscillator,
which corresponds to multiplication by $e^{-j(2\pi f_{c}t+\phi(t))}$ [59], the
equivalent complex baseband signal can be written as [44]

$\displaystyle y(t)$ $\displaystyle=\alpha\,s(t-\tau(t))e^{j\left[2\pi
f_{c}(t-\tau(t))+\phi(t-\tau(t))\right]}e^{-j\left[2\pi
f_{c}t+\phi(t)\right]}$ $\displaystyle=\alpha\,s(t-\tau(t))e^{-j2\pi
f_{c}\tau}e^{j2\pi f_{c}\nu
t}e^{j\left[\phi(t-\tau(t))-\phi(t)\right]}\leavevmode\nobreak\ .$ (5)

Focusing primarily on vehicular JRC scenarios, we assume that the Doppler
shifts satisfy $\lvert\nu\rvert\ll 1/N$ [38, 14, 20], where for typical
vehicular OFDM JRC systems, $\nu\ll 10^{-6}$ (corresponding to $v\ll
540\,\rm{km/h}$), while $N$ is on the order of $10^{3}$. This allows us to
approximate the PN term in (5) as $\phi(t-\tau(t))\approx\phi(t-\tau)$. In
addition, the time-bandwidth product $BMT_{\rm{sym}}$ is small enough to
justify (together with $\lvert\nu\rvert\ll 1/N$) the narrowband approximation
$s(t-\tau(t))\approx s(t-\tau)$ [51]. Under this setting, the received signal
in (5) becomes

$\displaystyle y(t)$ $\displaystyle=\alpha\,s(t-\tau)e^{-j2\pi
f_{c}\tau}e^{j2\pi f_{c}\nu t}w(t,\tau)\leavevmode\nobreak\ ,$ (6)

where the multiplicative PN process is represented by

$\displaystyle w(t,\tau)\triangleq
e^{j\left[\phi(t-\tau)-\phi(t)\right]}\leavevmode\nobreak\ .$ (7)

The statistical properties of the PN process $\phi(t-\tau)-\phi(t)$ in (7)
will be derived in Sec. III.

###  II-C Fast-Time/Slow-Time Representation with Phase Noise

For the ${m^{\rm{th}}}$ symbol, we remove the CP and sample $y(t)$ in (6) at
$t=mT_{\rm{sym}}+T_{\rm{cp}}+\ell T/N$ for $\ell=0,\ldots,N-1$. Making the
standard OFDM radar assumptions $T_{\rm{cp}}\geq\tau$ [38, 39, 1] (CP duration
is set to be longer than the round-trip delay of the furthermost target) and
$f_{c}T\nu\ll 1$ [60, 61, 62, 15] (Doppler shift $f_{c}\nu$ is small compared
to subcarrier spacing $\Delta f$), and ignoring constant phase terms, the
received signal for the ${m^{\rm{th}}}$ symbol can be written as [14, 20]

$\displaystyle y_{\ell,m}$ $\displaystyle=\alpha\,e^{j2\pi
f_{c}mT_{\rm{sym}}\nu}w_{\ell,m}(\tau)$ (8) $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \times\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x_{n,m}\,e^{j2\pi
n\frac{\ell}{N}}e^{-j2\pi n\Delta f\tau}\leavevmode\nobreak\ ,$

where $w_{\ell,m}(\tau)$ is the sampled version of the PN term $w(t,\tau)$ in
(6) for $t=mT_{\rm{sym}}+T_{\rm{cp}}+\ell T/N$. Let

$\displaystyle\mathbf{b}(\tau)$
$\displaystyle\triangleq{\left[1,e^{-j2\pi\Delta
f\tau},\ldots,e^{-j2\pi(N-1)\Delta f\tau}\right]}^{T}\leavevmode\nobreak\ ,$
(9) $\displaystyle\mathbf{c}(\nu)$ $\displaystyle\triangleq{\left[1,e^{-j2\pi
f_{c}T_{\rm{sym}}\nu},\ldots,e^{-j2\pi
f_{c}(M-1)T_{\rm{sym}}\nu}\right]}^{T}\leavevmode\nobreak\ ,$ (10)

represent the frequency-domain and temporal (slow-time) steering vectors,
respectively.

Aggregating the observations in (8) over fast-time $\ell$ and slow-time $m$,
and taking into account the presence of additive sensor noise, the fast-
time/slow-time observation matrix in the presence of PN is obtained as [20]

$\displaystyle\mathbf{Y}=\alpha\,\mathbf{W}\odot\mathbf{F}_{N}^{H}\Big{(}\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\Big{)}+\mathbf{Z}\leavevmode\nobreak\
,$ (11)

where $\mathbf{W}\in\mathbb{C}^{N\times M}$ with
$\left[\mathbf{W}\right]_{\ell,m}\triangleq w_{\ell,m}(\tau)$ is the
multiplicative PN matrix333For notational convenience, we drop the dependence
of $\mathbf{W}$ on $\tau$. consisting of fast-time/slow-time samples from the
PN process in (7), $\mathbf{F}_{N}\in\mathbb{C}^{N\times N}$ is the unitary
DFT matrix with
$\left[\mathbf{F}_{N}\right]_{\ell,n}=\frac{1}{\sqrt{N}}e^{-j2\pi
n\frac{\ell}{N}}$, $\mathbf{X}\in\mathbb{C}^{N\times M}$ contains the complex
data symbols with $\left[\mathbf{X}\right]_{n,m}\triangleq x_{n,m}$,
$\mathbf{Y}\in\mathbb{C}^{N\times M}$ with
$\left[\mathbf{Y}\right]_{\ell,m}\triangleq y_{\ell,m}$, and
$\mathbf{Z}\in\mathbb{C}^{N\times M}$ is additive white Gaussian noise (AWGN)
with
${\rm{vec}}\left(\mathbf{Z}\right)\sim{\mathcal{CN}}({{\bm{0}}},\allowbreak
2\sigma^{2}{\bm{\mathrm{I}}})$. As observed from (11), the PN component
$\mathbf{W}$ introduces intercarrier interference (ICI) in OFDM radar [51, 14,
20] (similar to its effect in OFDM communications [25, 19, 16, 21]) and might
severely degrade the performance of delay-Doppler estimation.

###  II-D Special Case: Ideal Oscillator

To relate the derived signal model in (11) to the commonly used ones in the
literature, we investigate the special case of an ideal oscillator where the
PN process is not present, i.e., $\phi(t)=0\,,\forall t$, which yields
$w(t,\tau)=1\,,\forall t,\tau$, and $\mathbf{W}$ becomes an all-ones matrix,
i.e., $\mathbf{W}={{\bm{1}}}_{N\times M}$. In this case, (11) reverts to the
PN-free model

$\displaystyle\mathbf{Y}^{\rm{PN-
free}}=\alpha\,\mathbf{F}_{N}^{H}\Big{(}\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\Big{)}+\mathbf{Z}\leavevmode\nobreak\
.$ (12)

Following the traditional processing chain for OFDM radar receivers [40, 12,
39], we take the DFT of the columns of $\mathbf{Y}^{\rm{PN-free}}$ in (12) to
switch from fast-time/slow-time to frequency/slow-time domain and obtain the
standard OFDM radar observations [15, 61, 60, 12, 40, 39]:

$\displaystyle\widetilde{\mathbf{Y}}=\mathbf{F}_{N}\mathbf{Y}^{\rm{PN-
free}}=\alpha\,\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)+\mathbf{F}_{N}\mathbf{Z}\leavevmode\nobreak\
,$ (13)

where
${\rm{vec}}\left(\mathbf{F}_{N}\mathbf{Z}\right)\sim{\mathcal{CN}}({{\bm{0}}},\allowbreak
2\sigma^{2}{\bm{\mathrm{I}}})$. Clearly, (13) does not involve any ICI effect,
and is therefore amenable to conventional delay-Doppler estimation algorithms
(after removing the effect of $\mathbf{X}$), such as 2-D DFT over time and
frequency domains [12, 39, 40] and super-resolution methods [60, 61].

###  II-E Problem Statement for OFDM Radar under Phase Noise

Given the transmit data symbols444Being co-located on a shared platform with
the JRC transmitter, the radar receiver has the knowledge of transmit data
symbols $\mathbf{X}$ [12, 14, 20]. $\mathbf{X}$ and the fast-time/slow-time
observations $\mathbf{Y}$ in (11), the problem of interest for OFDM radar
sensing in the presence of PN is to estimate the target parameters $\alpha$,
$\tau$ and $\nu$, which inherently involves estimating the PN matrix
$\mathbf{W}$ and compensating for its effect on $\mathbf{Y}$. To tackle this
problem, we first derive the statistical properties of the PN process in Sec.
III, which will then be utilized in Sec. IV to propose a novel algorithm for
joint estimation of delay, Doppler and PN. In Sec. V, we take a step further
by exploiting PN as something beneficial for radar sensing.

##  III Phase Noise Statistics

This section provides a statistical characterization of the PN process
$\phi(t-\tau)-\phi(t)$ in (7), gives expressions of the PN variance for
different types of oscillators and derives the structure of the PN covariance
matrix.

###  III-A Statistics of Differential Phase Noise Process

Let $\phi(t)$ be a zero-mean Gaussian random process with variance
$\sigma^{2}_{\phi}(t)$ [49, 30], i.e.,

$\phi(t)\sim{\mathcal{N}}(0,\sigma^{2}_{\phi}(t))\leavevmode\nobreak\ ,$ (14)

where the form of $\sigma^{2}_{\phi}(t)$ depends on the type of oscillator.
Let us define the differential PN (DPN) process [63] (also called self-
referenced PN, increment PN process [48], or PN variation [49]) as

$\xi(t,\tau)\triangleq\phi(t)-\phi(t-\tau)\leavevmode\nobreak\ .$ (15)

Since the DPN process555While $\xi(t,\tau)$ is called the DPN here to
distinguish it from the actual PN process $\phi(t)$, we will mostly refer to
$\xi(t,\tau)$ as PN in the remainder of the text for ease of exposition. is
stationary, its statistics depend only on the increment value (target delay)
$\tau$ [48, Sec. IV]. Hence, the DPN process can be statistically
characterized as [49]

$\xi(t,\tau)\sim{\mathcal{N}}(0,\sigma^{2}_{\xi}(\tau))\leavevmode\nobreak\ ,$
(16)

where $\sigma^{2}_{\xi}(\tau)$ is the delay-dependent variance of
$\xi(t,\tau)$. The following lemma provides the second-order statistics of
$\xi(t,\tau)$.

######  Lemma 1.

The correlation function of the DPN process $\xi(t,\tau)$ in (15) is given by

$\displaystyle R_{\xi\xi}(\Delta t,\tau)$
$\displaystyle=R_{\xi\xi}(t_{1},t_{2},\tau)\triangleq\operatorname*{\mathbb{E}}\left[\xi(t_{1},\tau)\xi(t_{2},\tau)\right]$
(17) $\displaystyle=\frac{\sigma^{2}_{\xi}(\tau+\Delta
t)+\sigma^{2}_{\xi}(\tau-\Delta t)}{2}-\sigma^{2}_{\xi}(\Delta
t)\leavevmode\nobreak\ ,$

where $\Delta t\triangleq t_{1}-t_{2}$ is the time difference.

######  Proof.

Please see Appendix A.

###  III-B PN Variance for Different Oscillator Types

Lemma 1 allows us to compute the correlation of $\xi(t,\tau)$ in terms of its
delay-dependent variance function $\sigma^{2}_{\xi}(\tau)$. We now provide
expressions for $\sigma^{2}_{\xi}(\tau)$ by distinguishing between two
realizations of an oscillator, namely, free-running oscillator (FRO) and
phase-locked loop (PLL) synthesizer [25, 48].

####  III-B1 Free-Running Oscillators (FROs)

For FROs, the variance of the DPN $\xi(t,\tau)$ in (16) is given by [48, Sec.
V], [49, Sec. V, VI], [25, Sec. III-A]

$\displaystyle\sigma^{2}_{\xi}(\tau)$ $\displaystyle=4\pi
f_{\rm{3dB}}\lvert\tau\rvert\leavevmode\nobreak\ ,$ (18)

where $f_{\rm{3dB}}$ is the $3\,\rm{dB}$ bandwidth of the Lorentzian
oscillator spectrum.

####  III-B2 Phase-Locked Loop (PLL) Synthesizers

For PLL architectures, the delay-dependent variance in (16) can be expressed
as [48, Sec. VII-A], [49, Sec. III]

$\displaystyle\sigma^{2}_{\xi}(\tau)$
$\displaystyle=\frac{2f_{\rm{3dB}}}{f_{\rm{loop}}}\left(1-e^{-2\pi
f_{\rm{loop}}\lvert\tau\rvert}\right)\leavevmode\nobreak\ ,$ (19)

where $f_{\rm{loop}}$ denotes the loop bandwidth of PLL.

We note that PLL degenerates to FRO with decreasing $f_{\rm{loop}}$, i.e.,
(19) converges to (18) as $f_{\rm{loop}}\to 0$. With the expressions (18) and
(19), a complete statistical characterization of $\xi(t,\tau)$ can be obtained
via (16) and (17) for FROs and PLLs, which, in turn, yields the covariance
matrix of the fast-time/slow-time PN samples in $\mathbf{W}$ in (11). As an
example, Fig. 2 plots the covariance of $\xi(t,\tau)$ at different target
delays $\tau$ for both PLLs and FROs. The figure illustrates the delay-
dependency of the PN statistics in OFDM radar sensing (as opposed to OFDM
communications, e.g., [46, 47, 19]) and provides insights into the effect of
PLL based control on the correlation behavior of PN.

Figure 2: Covariance of the DPN process $\xi(t,\tau)$ in (15), computed via
the expression in (17), for PLL and FRO architectures at two different target
ranges $R=50\,\rm{m}$ and $R=150\,\rm{m}$, where the $3\,\rm{dB}$ bandwidth is
$f_{\rm{3dB}}=50\,\rm{kHz}$ and the PLL loop bandwidth is
$f_{\rm{loop}}=500\,\rm{kHz}$. Contrary to OFDM communications, the PN
statistics in OFDM radar depend on target delay, as also illustrated in Fig.
1.

###  III-C Delay-Dependent PN Covariance Matrix

Let
$\mathbf{w}\triangleq{\rm{vec}}\left(\mathbf{W}\right)\in\mathbb{C}^{NM\times
1}$ and $\bm{\xi}\in\mathbb{R}^{NM\times 1}$ be the sampled version of the PN
process in (15) over the entire OFDM frame, i.e., $\mathbf{w}=e^{-j\bm{\xi}}$
from (7). Then, $\bm{\xi}$ can be statistically characterized as666The
derivation of PN statistics can be straightforwardly extended to the multi-
target case by computing cross-correlation of PN vectors associated with
different targets at different delays, following the arguments in Lemma 1.

$\displaystyle\bm{\xi}\sim{\mathcal{N}}({{\bm{0}}},\mathbf{R}(\tau))\leavevmode\nobreak\
,$ (20)

where $\mathbf{R}(\tau)\in\mathbb{R}^{NM\times NM}$ is the delay-dependent
positive definite covariance matrix of $\bm{\xi}$. Using (16) and (17), the
${(i_{1},i_{2})^{\rm{th}}}$ entry of $\mathbf{R}(\tau)$ can be written as

$\displaystyle\left[\mathbf{R}(\tau)\right]_{i_{1},i_{2}}=R_{\xi\xi}(\Delta
t_{i_{1}i_{2}},\tau)\leavevmode\nobreak\ ,$ (21)

where

$\displaystyle\Delta t_{i_{1}i_{2}}$
$\displaystyle\triangleq(i_{1}-i_{2})T_{\rm{s}}+(m_{1}-m_{2})T_{\rm{cp}}$ (22)

for $(i_{1},i_{2})=(n_{1}+m_{1}N,n_{2}+m_{2}N)$, with $0\leq n_{1},n_{2}\leq
N-1$ and $0\leq m_{1},m_{2}\leq M-1$ denoting the fast-time and slow-time
sample indices, respectively, and $T_{\rm{s}}=T/N$ the sampling interval. We
note from (22) that due to CP removal, $\Delta t_{i_{1}i_{2}}$ depends not
only on $i_{1}-i_{2}$, but also on the difference between symbol (slow-time)
indices for PN samples belonging to different symbols.

From (21) and (22), it is straightforward to see that $\mathbf{R}(\tau)$ is a
symmetric Toeplitz-block Toeplitz matrix [64] consisting of $M\times M$ blocks
of size $N\times N$:

$\displaystyle\mathbf{R}(\tau)=\begin{bmatrix}\mathbf{R}_{0}(\tau)&\mathbf{R}_{1}(\tau)&\ldots&\mathbf{R}_{M-1}(\tau)\\\
\mathbf{R}_{1}^{T}(\tau)&\ddots&\ddots&\\\
&\ddots&\ddots&\mathbf{R}_{1}(\tau)\\\
\mathbf{R}_{M-1}^{T}(\tau)&&\mathbf{R}_{1}^{T}(\tau)&\mathbf{R}_{0}(\tau)\end{bmatrix}\leavevmode\nobreak\
,$ (23)

where the ${m^{\rm{th}}}$ Toeplitz block
$\mathbf{R}_{m}(\tau)\in\mathbb{R}^{N\times N}$ is given by

$\displaystyle[\mathbf{R}_{m}(\tau)]_{n_{1},n_{2}}=R_{\xi\xi}(\Delta
t^{(m)}_{n_{1}n_{2}},\tau)$ (24)

with $\Delta
t^{(m)}_{n_{1}n_{2}}\triangleq(n_{1}-n_{2})T_{\rm{s}}-mT_{\rm{sym}}$. While
PLLs exhibit the generic Toeplitz-block Toeplitz structure in (23), FROs lead
to a more special covariance structure, as pointed out in the following lemma.

######  Lemma 2.

The PN covariance matrix $\mathbf{R}(\tau)$ in (23) is block-diagonal for FROs
with $\tau\leq T_{\rm{cp}}$, i.e.,

$\displaystyle\mathbf{R}(\tau)={\rm{blkdiag}}\big{(}\mathbf{R}_{0}(\tau),\ldots,\mathbf{R}_{0}(\tau)\big{)}\leavevmode\nobreak\
.$ (25)

######  Proof.

Please see Appendix B.

An intuitive interpretation of Lemma 2 can be provided as follows. Since PN
samples in different symbols are separated in time by at least $T_{\rm{cp}}$,
the PN process in (15) with $\tau\leq T_{\rm{cp}}$ becomes uncorrelated from
one symbol to another. In other words, the time intervals for PN accumulation
[48] corresponding to $\xi(t,\tau)$ and $\xi(t+T_{\rm{cp}},\tau)$ are non-
overlapping, leading to uncorrelated PN in the absence of a control loop. This
result is also corroborated by Fig. 2, where the correlation in case of FRO is
zero for $\Delta t\geq\tau$, while the correlation for PLL can have non-zero
values for $\Delta t\geq T_{\rm{cp}}$.

##  IV Proposed Algorithm for Delay-Doppler Estimation under Phase Noise

In this section, using the statistical characterization of PN derived in (20)
and (23), we formulate the sensing problem stated in Sec. II-E using a hybrid
ML/MAP estimation approach and propose a novel iterated small angle
approximation (ISAA) algorithm to jointly estimate delay, Doppler and PN.

###  IV-A Hybrid ML/MAP Estimator

To derive the hybrid ML/MAP estimator of delay, Doppler and PN, we first
rewrite the observations in (11) as

$\displaystyle\mathbf{y}=\alpha\,\mathbf{\Xi}\mathbf{q}(\tau,\nu)+\mathbf{z}\leavevmode\nobreak\
,$ (26)

where
$\mathbf{y}\triangleq{\rm{vec}}\left(\mathbf{Y}\right)\in\mathbb{C}^{NM\times
1}$,
$\mathbf{z}\triangleq{\rm{vec}}\left(\mathbf{Z}\right)\in\mathbb{C}^{NM\times
1}$,

$\displaystyle\mathbf{\Xi}$
$\displaystyle\triangleq{\rm{diag}}\left(e^{-j\bm{\xi}}\right)\in\mathbb{C}^{NM\times
NM}\leavevmode\nobreak\ ,$ (27) $\displaystyle\mathbf{q}(\tau,\nu)$
$\displaystyle\triangleq{\rm{vec}}\left(\mathbf{F}_{N}^{H}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{[}}\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{]}}\right)\in\mathbb{C}^{NM\times
1}\leavevmode\nobreak\ .$ (28)

Our goal herein is to estimate from (26) the unknown parameter vector
${\bm{\eta}}=\left[\tau,\nu,\alpha,\bm{\xi}^{T}\right]^{T}$, consisting of
both random ($\bm{\xi}$) and deterministic ($\tau,\nu,\alpha$)
parameters777Estimation of ${\bm{\eta}}$ should be performed for each new OFDM
frame due to the time-varying nature of the PN and possible movement of the
target in the delay-Doppler plane. Although the PN $\bm{\xi}$ has time-
invariant statistics specified in (20), it will have a different (time-
varying) realization in each frame and thus needs to be re-estimated by re-
executing the proposed method in Sec. IV as new observations in the form of
(26) arrive.. Then, the hybrid ML/MAP estimator of ${\bm{\eta}}$ can be
written as [53]

$\displaystyle({\widehat{\tau}},{\widehat{\nu}},\widehat{\alpha},\widehat{\bm{\xi}})=\arg\max_{\tau,\nu,\alpha,\bm{\xi}}f_{\mathbf{y},\bm{\xi}}(\mathbf{y},\bm{\xi};\tau,\nu,\alpha),$
(29)

where $f_{\mathbf{y},\bm{\xi}}(\mathbf{y},\bm{\xi};\tau,\nu,\alpha)$ is the
joint PDF of $\mathbf{y}$ and $\bm{\xi}$ as a function of the deterministic
parameters $\tau$, $\nu$ and $\alpha$

$\displaystyle
f_{\mathbf{y},\bm{\xi}}(\mathbf{y},\bm{\xi};\tau,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu},\alpha)=f_{\mathbf{y}\lvert\bm{\xi}}(\mathbf{y}\lvert\bm{\xi};\tau,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu},\alpha)f_{\bm{\xi}}(\bm{\xi};\tau)\leavevmode\nobreak\
,$ (30)

$f_{\mathbf{y}\lvert\bm{\xi}}(\mathbf{y}\lvert\bm{\xi};\tau,\nu,\alpha)$ is
the conditional PDF of $\mathbf{y}$ given $\bm{\xi}$, and
$f_{\bm{\xi}}(\bm{\xi};\tau)$ is the a-priori PDF of $\bm{\xi}$. It follows
from (20) and (26) that

$\displaystyle
f_{\mathbf{y}\lvert\bm{\xi}}(\mathbf{y}\lvert\bm{\xi};\tau,\nu,\alpha)$ (31)
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{NM}}\exp\left\\{-\frac{\left\lVert\mathbf{y}-\alpha\,\mathbf{\Xi}\mathbf{q}(\tau,\nu)\right\rVert^{2}}{2\sigma^{2}}\right\\}\leavevmode\nobreak\
,$ $\displaystyle f_{\bm{\xi}}(\bm{\xi};\tau)$ (32)
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{NM}{{{\rm{det}}\left(\mathbf{R}(\tau)\right)}}}}\exp\left\\{-\frac{\bm{\xi}^{T}\mathbf{R}(\tau)^{-1}\bm{\xi}}{2}\right\\}\leavevmode\nobreak\
.$

Plugging (30)–(32) into (29) yields

$\displaystyle({\widehat{\tau}},{\widehat{\nu}},\widehat{\alpha},\widehat{\bm{\xi}})$
$\displaystyle=\arg\min_{\tau,\nu,\alpha,\bm{\xi}}\Bigg{\\{}\frac{\left\lVert\mathbf{y}-\alpha\,\mathbf{\Xi}\mathbf{q}(\tau,\nu)\right\rVert^{2}}{2\sigma^{2}}$
(33) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\frac{\bm{\xi}^{T}\mathbf{R}(\tau)^{-1}\bm{\xi}+\log\det\mathbf{R}(\tau)}{2}\Bigg{\\}}\leavevmode\nobreak\
.$

For given $\tau$, $\nu$ and $\bm{\xi}$, the optimal estimate of $\alpha$ in
(33) is given by

$\displaystyle\widehat{\alpha}=\frac{\mathbf{q}^{H}(\tau,\nu)\mathbf{\Xi}^{H}\mathbf{y}}{\mathbf{q}^{H}(\tau,\nu)\mathbf{\Xi}^{H}\mathbf{\Xi}\mathbf{q}(\tau,\nu)}=\frac{\mathbf{q}^{H}(\tau,\nu)\mathbf{\Xi}^{H}\mathbf{y}}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\leavevmode\nobreak\
,$ (34)

where the last equality follows from (9), (10) and (28):

$\displaystyle\left\lVert\mathbf{q}(\tau,\nu)\right\rVert^{2}$
$\displaystyle=\left\lVert\mathbf{F}_{N}^{H}\Big{(}\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\Big{)}\right\rVert_{F}^{2}$
(35)
$\displaystyle=\left\lVert\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\right\rVert_{F}^{2}=\left\lVert\mathbf{X}\right\rVert_{F}^{2}\leavevmode\nobreak\
.$

Inserting (34) into (33), the hybrid ML/MAP problem becomes

$\displaystyle({\widehat{\tau}},{\widehat{\nu}},\widehat{\bm{\xi}})$
$\displaystyle=\arg\min_{\tau,\nu,\bm{\xi}}\leavevmode\nobreak\
\frac{\mathbf{y}^{H}\bm{\Pi}^{\perp}_{\mathbf{\Xi}\mathbf{q}(\tau,\nu)}\mathbf{y}}{\sigma^{2}}+\bm{\xi}^{T}\mathbf{R}(\tau)^{-1}\bm{\xi}$
(36) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\log\det\mathbf{R}(\tau)\leavevmode\nobreak\ ,$

where

$\displaystyle\bm{\Pi}^{\perp}_{\mathbf{\Xi}\mathbf{q}(\tau,\nu)}$
$\displaystyle={\bm{\mathrm{I}}}-\frac{\mathbf{\Xi}\mathbf{q}(\tau,\nu)\mathbf{q}^{H}(\tau,\nu)\mathbf{\Xi}^{H}}{\left\lVert\mathbf{q}(\tau,\nu)\right\rVert^{2}}$
(37)
$\displaystyle=\mathbf{\Xi}\mathbf{\Xi}^{H}-\frac{\mathbf{\Xi}\mathbf{q}(\tau,\nu)\mathbf{q}^{H}(\tau,\nu)\mathbf{\Xi}^{H}}{\left\lVert\mathbf{q}(\tau,\nu)\right\rVert^{2}}$
$\displaystyle=\mathbf{\Xi}\left({\bm{\mathrm{I}}}-\frac{\mathbf{q}(\tau,\nu)\mathbf{q}^{H}(\tau,\nu)}{\left\lVert\mathbf{q}(\tau,\nu)\right\rVert^{2}}\right)\mathbf{\Xi}^{H}$
$\displaystyle=\mathbf{\Xi}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}\mathbf{\Xi}^{H}\leavevmode\nobreak\
,$

with $\mathbf{\Xi}\mathbf{\Xi}^{H}={\bm{\mathrm{I}}}$ resulting from (27).
Substituting (37) into (36), we obtain

$\displaystyle({\widehat{\tau}},{\widehat{\nu}},\widehat{\bm{\xi}})=\arg\min_{\tau,\nu,\bm{\xi}}\leavevmode\nobreak\
\mathcal{L}(\tau,\nu,\bm{\xi})\leavevmode\nobreak\ ,$ (38)

where

$\displaystyle\mathcal{L}(\tau,\nu,\bm{\xi})$
$\displaystyle=\frac{\mathbf{y}^{H}\mathbf{\Xi}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}\mathbf{\Xi}^{H}\mathbf{y}}{\sigma^{2}}+\bm{\xi}^{T}\mathbf{R}(\tau)^{-1}\bm{\xi}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\log\det\mathbf{R}(\tau)\leavevmode\nobreak\ .$ (39)

###  IV-B Iterated Small Angle Approximation

The hybrid ML/MAP optimization problem derived in (38) seems quite challenging
to solve, primarily due to highly non-linear behavior of the high-dimensional
PN vector $\bm{\xi}$ in the objective (39) (see (27)), leading to many local
optima. A possible remedy to overcome such non-linearity is to employ small
angle approximation (SAA) $e^{j\theta}\approx 1+j\theta$ [29] for small
$\theta$. However, while this approach can work well for PN estimation in
communications, it may lead to large errors in our sensing problem of
interest. This is especially the case for distant targets since the PN
variance increases with target delay, as seen from (18) and (19), invalidating
the assumption of small $\theta$.

To circumvent the non-linearity of PN in (39) and deal with large PN
variances, we propose an iterated small angle approximation (ISAA) approach
that invokes SAA around a given estimate of PN at each iteration, starting
from an all-zeros estimate $\widehat{\bm{\xi}}={{\bm{0}}}$. More specifically,
suppose we have an estimate of PN vector at the ${i^{\rm{th}}}$ iteration,
denoted by $\widehat{\bm{\xi}}^{(i)}\in\mathbb{R}^{NM\times 1}$, and we wish
to approximate the exponential PN term $e^{-j\bm{\xi}}$ in (39) around
$\widehat{\bm{\xi}}^{(i)}$. To this end, we define the residual PN

$\displaystyle\bm{\xi}_{\Delta}=\bm{\xi}-\widehat{\bm{\xi}}^{(i)}$ (40)

for a given PN estimate $\widehat{\bm{\xi}}^{(i)}$. Re-writing (39) as a
function of $\bm{\xi}_{\Delta}$, we have

$\displaystyle\mathcal{L}(\tau,\nu,\bm{\xi}_{\Delta}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\widehat{\bm{\xi}}^{(i)}})=\frac{(\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}})^{T}{\rm{diag}}\left(\mathbf{y}\right)^{H}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}{\rm{diag}}\left(\mathbf{y}\right)(\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}})^{\ast}}{\sigma^{2}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+(\bm{\xi}_{\Delta}+\widehat{\bm{\xi}}^{(i)})^{T}\mathbf{R}(\tau)^{-1}(\bm{\xi}_{\Delta}+\widehat{\bm{\xi}}^{(i)})+\log\det\mathbf{R}(\tau)\leavevmode\nobreak\
,$ (41)

where
$\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}=e^{-j(\bm{\xi}_{\Delta}+\widehat{\bm{\xi}}^{(i)})}\in\mathbb{C}^{NM\times
1}$. In order to solve (38) for $\bm{\xi}$ in an iterative fashion, we propose
to solve the following minimization problem at the ${i^{\rm{th}}}$ iteration
for the residual PN $\bm{\xi}_{\Delta}$, given a PN estimate
$\widehat{\bm{\xi}}^{(i)}$ at hand:

$\displaystyle({\widehat{\tau}},{\widehat{\nu}},\widehat{\bm{\xi}}_{\Delta})=\arg\min_{\tau,\nu,\bm{\xi}_{\Delta}}\leavevmode\nobreak\
\mathcal{L}(\tau,\nu,\bm{\xi}_{\Delta}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\widehat{\bm{\xi}}^{(i)}})\leavevmode\nobreak\
.$ (42)

To tackle (42), we invoke SAA for
$\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}$
in (41) around $\bm{\xi}_{\Delta}$ to obtain

$\displaystyle\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}=e^{-j\widehat{\bm{\xi}}^{(i)}}\odot
e^{-j\bm{\xi}_{\Delta}}\approx
e^{-j\widehat{\bm{\xi}}^{(i)}}\odot({{\bm{1}}}-j\bm{\xi}_{\Delta})\leavevmode\nobreak\
.$ (43)

Plugging (43) into the first term in (41), we obtain the approximation

$\displaystyle(\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}})^{T}{\rm{diag}}\left(\mathbf{y}\right)^{H}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}{\rm{diag}}\left(\mathbf{y}\right)(\bm{\upsilon}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}})^{\ast}$
$\displaystyle\approx\left[e^{-j\widehat{\bm{\xi}}^{(i)}}\odot({{\bm{1}}}-j\bm{\xi}_{\Delta})\right]^{T}{\rm{diag}}\left(\mathbf{y}\right)^{H}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}{\rm{diag}}\left(\mathbf{y}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\times\left[e^{j\widehat{\bm{\xi}}^{(i)}}\odot({{\bm{1}}}+j\bm{\xi}_{\Delta})\right]$
$\displaystyle=({{\bm{1}}}-j\bm{\xi}_{\Delta})^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)({{\bm{1}}}+j\bm{\xi}_{\Delta})\leavevmode\nobreak\
,$ (44)

where
$\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\in\mathbb{C}^{NM\times
NM}$ is defined as888Note that (44) is guaranteed to be real since
$\mathbf{\Gamma}(\tau,\nu)$ in (45) is Hermitian.

$\displaystyle\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)$
$\displaystyle\triangleq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{[}}{\rm{diag}}\left(\mathbf{y}\right)^{H}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}{\rm{diag}}\left(\mathbf{y}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{]}}\odot{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{[}}e^{-j\widehat{\bm{\xi}}^{(i)}}(e^{j\widehat{\bm{\xi}}^{(i)}})^{T}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{]}}\leavevmode\nobreak\
.$ (45)

Now, substituting the approximation (44) into (41), we have

$\displaystyle\mathcal{L}(\tau,\nu,\bm{\xi}_{\Delta}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\widehat{\bm{\xi}}^{(i)}})$
(46)
$\displaystyle\approx\frac{1}{\sigma^{2}}({{\bm{1}}}-j\bm{\xi}_{\Delta})^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)({{\bm{1}}}+j\bm{\xi}_{\Delta})$
(47) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+(\widehat{\bm{\xi}}^{(i)})^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}+\bm{\xi}_{\Delta}^{T}\mathbf{R}(\tau)^{-1}\bm{\xi}_{\Delta}+2\bm{\xi}_{\Delta}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\log\det\mathbf{R}(\tau)$
$\displaystyle=\frac{1}{\sigma^{2}}\Big{[}{{\bm{1}}}^{T}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+\bm{\xi}_{\Delta}^{T}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\bm{\xi}_{\Delta}$
(48) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+2\bm{\xi}_{\Delta}^{T}\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu\right\\}{{\bm{1}}}\Big{]}+\bm{\xi}_{\Delta}^{T}\mathbf{R}(\tau)^{-1}\bm{\xi}_{\Delta}+2\bm{\xi}_{\Delta}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+(\widehat{\bm{\xi}}^{(i)})^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}+\log\det\mathbf{R}(\tau)$
$\displaystyle=\bm{\xi}_{\Delta}^{T}\left(\frac{1}{\sigma^{2}}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}+\mathbf{R}(\tau)^{-1}\right)\bm{\xi}_{\Delta}$
(49) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+2\bm{\xi}_{\Delta}^{T}\left(\frac{1}{\sigma^{2}}\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{1}{\sigma^{2}}{{\bm{1}}}^{T}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+(\widehat{\bm{\xi}}^{(i)})^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}+\log\det\mathbf{R}(\tau)\leavevmode\nobreak\
,$

where (48) is due to
$\bm{\xi}_{\Delta}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}$
being a real vector and $\mathbf{\Gamma}(\tau,\nu)$ being a Hermitian matrix
(please see Appendix C for details).

Observing that (49) is quadratic in $\bm{\xi}_{\Delta}$, the optimal estimate
of $\bm{\xi}_{\Delta}$ that minimizes the approximated version of
$\mathcal{L}(\tau,\nu,\bm{\xi}_{\Delta}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\widehat{\bm{\xi}}^{(i)}})$
in (49) can be written in closed-form for a given delay-Doppler pair
$(\tau,\nu)$ as follows:

$\displaystyle\widehat{\bm{\xi}}_{\Delta}(\tau,\nu)$
$\displaystyle=-\mathbf{R}(\tau)\Big{(}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)+\sigma^{2}{\bm{\mathrm{I}}}\Big{)}^{-1}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\times\left(\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+\sigma^{2}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}\right)\leavevmode\nobreak\
.$ (50)

Finally, using the residual estimate in (50) and the definition in (40), the
PN estimate can be updated as

$\displaystyle\widehat{\bm{\xi}}^{(i+1)}=\widehat{\bm{\xi}}^{(i)}+\widehat{\bm{\xi}}_{\Delta}(\tau,\nu)\leavevmode\nobreak\
.$ (51)

###  IV-C Alternating Optimization to Solve (38)

The iterative procedure developed in (50) and (51) for updating the PN
estimate $\widehat{\bm{\xi}}$ as a function of delay and Doppler motivates an
alternating optimization method to solve the original hybrid ML/MAP
optimization problem (38). Hence, we propose to estimate delay, Doppler and PN
in (38) using an iterative refinement approach that alternates between PN
estimation and delay-Doppler estimation as follows:

  * •

Update $\widehat{\bm{\xi}}$ for Fixed $(\tau,\nu)$: For a given delay-Doppler
pair $(\tau,\nu)$, we compute the residual PN via (50) and update
$\widehat{\bm{\xi}}$ via (51).

  * •

Update $(\tau,\nu)$ for Fixed $\widehat{\bm{\xi}}$: For a given PN estimate
$\widehat{\bm{\xi}}$, we find the optimal delay-Doppler pair by solving (38):

$\displaystyle({\widehat{\tau}},{\widehat{\nu}})=\arg\min_{\tau,\nu}\leavevmode\nobreak\
\mathcal{L}(\tau,\nu,\widehat{\bm{\xi}})\leavevmode\nobreak\ ,$ (52)

which is equivalent to (please see Appendix D for details)

$\displaystyle({\widehat{\tau}},{\widehat{\nu}})$
$\displaystyle=\arg\max_{\tau,\nu}\leavevmode\nobreak\
\widetilde{\mathcal{L}}(\tau,\nu,\widehat{\bm{\xi}})\leavevmode\nobreak\ ,$
(53)

where

$\displaystyle\widetilde{\mathcal{L}}(\tau,\nu,\widehat{\bm{\xi}})$
$\displaystyle=\frac{\Big{\lvert}\mathbf{b}^{H}(\tau){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{[}}\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\big{(}\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\big{)}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Big{]}}\mathbf{c}(\nu)\Big{\rvert}^{2}}{\sigma^{2}\left\lVert\mathbf{X}\right\rVert_{F}^{2}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-\widehat{\bm{\xi}}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}-\log\det\mathbf{R}(\tau)\leavevmode\nobreak\
,$ (54)

and
$\widehat{\mathbf{W}}\triangleq{\rm{reshape}}_{N,M}\left(e^{-j\widehat{\bm{\xi}}}\right)$.

The entire algorithm to solve (38) is summarized in Algorithm 1 and referred
to as MAP-ISAA999Algorithm 1 is agnostic to data symbols $\mathbf{X}$ in the
sense that it can work well with arbitrary $\mathbf{X}$ and imposes no
constraints on data symbols generated by the communications system.. Under
certain conditions, Algorithm 1 converges to a stationary point of (38)
(please see Appendix E for a detailed convergence analysis). Employing the
conjugate gradient (CG) method [29] to evaluate (50), the per-iteration
complexity of Algorithm 1 can be obtained as (please see Appendix F for
details)

$\displaystyle\mathcal{O}\big{(}MN\left(\log M+(IM_{0}+1)\log
N\right)\big{)}\leavevmode\nobreak\ ,$ (55)

where $I$ is the number of CG iterations and $M_{0}$ is the number of dominant
blocks of $\mathbf{R}(\tau)$ in (23), with $M_{0}=1$ for FROs due to Lemma 2
and $1\leq M_{0}\ll M$ for PLLs (typically, $M_{0}\leq 3$).

Algorithm 1 Joint Delay, Doppler and PN Estimation with MAP Criterion via
Iterated Small Angle Approximation (MAP-ISAA).

1:Input: Fast-time/slow-time observations $\mathbf{Y}$ in (11), convergence
thresholds $\epsilon_{\tau}$ and $\epsilon_{\nu}$ for delay and Doppler, and
maximum number of iterations $I_{\rm{max}}$.

2:Output: Estimates of delay, Doppler and PN
$\\{{\widehat{\tau}},{\widehat{\nu}},\widehat{\bm{\xi}}\\}$.

3:Initialization: Set $i=0$.

  1. 1.

Initialize the PN estimate to be the all-zeros vector, i.e.,
$\widehat{\bm{\xi}}^{(0)}={{\bm{0}}}$, in accordance with (20).

  2. 2.

Initialize the delay-Doppler pair to be the output of the standard 2-D FFT
method, i.e.,

$\displaystyle({\widehat{\tau}}^{(0)},{\widehat{\nu}}^{(0)})$
$\displaystyle=\arg\max_{\tau,\nu}\leavevmode\nobreak\
\Big{\lvert}\mathbf{b}^{H}(\tau)\left(\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\mathbf{Y}\right)\mathbf{c}(\nu)\Big{\rvert}^{2}\leavevmode\nobreak\
.$

4:Iterated Approximation Steps:

5:while $i<I_{\rm{max}}$

6: Update PN estimate via (50) and (51):
$\displaystyle\widehat{\bm{\xi}}^{(i+1)}=\widehat{\bm{\xi}}^{(i)}+\widehat{\bm{\xi}}_{\Delta}({\widehat{\tau}}^{(i)},{\widehat{\nu}}^{(i)})\leavevmode\nobreak\
.$

7: Update delay-Doppler estimate via (53) and (54):
$\displaystyle({\widehat{\tau}}^{(i+1)},{\widehat{\nu}}^{(i+1)})$
$\displaystyle=\arg\max_{\tau,\nu}\leavevmode\nobreak\
\widetilde{\mathcal{L}}(\tau,\nu,\widehat{\bm{\xi}}^{(i+1)})\leavevmode\nobreak\
.$

8: Set $i=i+1$.

9: if
$\lvert{\widehat{\tau}}^{(i)}-{\widehat{\tau}}^{(i-1)}\rvert\leq\epsilon_{\tau}$
and
$\lvert{\widehat{\nu}}^{(i)}-{\widehat{\nu}}^{(i-1)}\rvert\leq\epsilon_{\nu}$

10: break

11: end if

12:end while

##  V Extension to Range-Ambiguous Targets: From Mitigation to Exploitation

We have up to now focused on how to estimate and mitigate PN in the OFDM
sensing problem of interest. In this section, we extend the proposed
estimation framework to the case of range-ambiguous targets and devise a PN
exploitation approach to resolve range ambiguity.

###  V-A PN Exploitation: Concept Description

An important peculiarity of PN in OFDM monostatic sensing is that it provides
additional source of information on target delays that can be exploited to
improve delay estimation performance, whereas PN in OFDM communications always
degrades performance (e.g., [46, 47, 19, 16, 21]). More precisely, it is seen
from (9), (20) and (28) that both $\mathbf{b}(\tau)$ and $\mathbf{R}(\tau)$,
which are functions of the unknown delay $\tau$, have an impact on the
observation (26) and the hybrid ML/MAP cost function (39). Hence, in addition
to the standard frequency-domain phase rotations in $\mathbf{b}(\tau)$, the
covariance matrix $\mathbf{R}(\tau)$ of the PN vector conveys information on
target delay. As noticed from the elements of $\mathbf{R}(\tau)$ in (17)–(21),
there is no maximum unambiguous delay imposed by $\mathbf{R}(\tau)$ in
estimating $\tau$. This implies that we can detect the true ranges of the
targets with $\tau>T$ by exploiting the information in $\mathbf{R}(\tau)$.
However, for the ideal case without PN, we can only extract the delay
information from $\mathbf{b}(\tau)$ (e.g., [60, 12, 40, 39]), which leads to
range ambiguity for targets with $\tau>T$ due to the periodicity of the
complex exponential terms in (9).

###  V-B PN Exploitation: Proposed Algorithm

Let us consider a range-ambiguous target with true (unambiguous) delay
$\tau=\tau_{\rm{p}}+kT$ for some integer $k\geq 1$ denoting the ambiguity
index, where $0<\tau_{\rm{p}}\leq T$ is the principal delay of the target,
i.e., $\tau_{\rm{p}}=\tau\ (\mathrm{mod}\ T)$. Based on the observations in
Sec. V-A, we propose to resolve the range ambiguity using the PN statistics
represented by $\mathbf{R}(\tau)$. The key insight here is that the
observation-related (first) term in (39) assumes the same value for $\tau$ and
$\tau_{\rm{p}}$ due to inherent range ambiguity in $\mathbf{b}(\cdot)$, while
the statistics of the PN estimate $\widehat{\bm{\xi}}$ at the output of
Algorithm 1 would match only with $\mathbf{R}(\tau)$ (not with
$\mathbf{R}(\tau_{\rm{p}})$), enabling range ambiguity resolution101010Please
see Fig. 11 in Sec. VI-E for an illustration of this phenomenon.. Inspired by
this insight, we formulate a parametric covariance matrix reconstruction
problem111111The covariance matching problem formulated in (56) has a
theoretical justification based on the extended invariance principle (EXIP)
[65, 66] and the ML estimator of delay from $\widehat{\bm{\xi}}$. Please refer
to Appendix G for details. that minimizes the Frobenius distance between the
sample matrix $\widehat{\mathbf{R}}=\widehat{\bm{\xi}}\widehat{\bm{\xi}}^{T}$
and the parametric matrix $\mathbf{R}(\tau)$ in (23):

$\displaystyle{{\widehat{\tau}}^{\rm{true}}}=\arg\min_{\tau}\leavevmode\nobreak\
\big{\lVert}\mathbf{R}(\tau)-\widehat{\mathbf{R}}\big{\rVert}_{F}^{2}\leavevmode\nobreak\
,$ (56)

potentially yielding the unambiguous (true) delay estimate
${{\widehat{\tau}}^{\rm{true}}}$. Exploiting the Toeplitz-block Toeplitz
structure in (23), the problem (56) reduces to [67]

$\displaystyle{{\widehat{\tau}}^{\rm{true}}}=\arg\min_{\tau}\leavevmode\nobreak\
\big{\lVert}[\mathbf{R}(\tau)]_{0,:}-\widehat{\mathbf{r}}\big{\rVert}^{2}\leavevmode\nobreak\
,$ (57)

where $[\mathbf{R}(\tau)]_{0,:}$ is the first row of $\mathbf{R}(\tau)$ and

$\displaystyle\widehat{r}_{i}=\frac{1}{NM-i}\sum_{k=0}^{NM-i-1}\widehat{\xi}_{k}\,\widehat{\xi}_{k+i}\leavevmode\nobreak\
,$ (58)

with $\widehat{\mathbf{r}}=[\widehat{r}_{0},\ldots,\widehat{r}_{NM-1}]^{T}$
and $\widehat{\bm{\xi}}=[\widehat{\xi}_{0},\ldots,\widehat{\xi}_{NM-1}]^{T}$.

In Algorithm 2, we present the proposed algorithm for range ambiguity
resolution via PN exploitation. We restrict the search interval for range in
Algorithm 1 such that the resulting delay estimate ${\widehat{\tau}}$ is
ambiguous. As seen from Algorithm 2, the prior information extracted from the
PN estimates is exploited to find the range ambiguity index of the target
(i.e., ambiguity resolution), while the principal range estimate of the target
comes from Algorithm 1.

Algorithm 2 PN Exploitation to Resolve Range Ambiguity.

1:Input: Estimates of (ambiguous) delay and PN at the output of Algorithm 1
$\\{{\widehat{\tau}},\widehat{\bm{\xi}}\\}$.

2:Output: Unambiguous (true) estimate of delay
${{\widehat{\tau}}^{\rm{true}}}$.

  1. 1.

Compute an estimate
$\widehat{\mathbf{r}}=[\widehat{r}_{0},\ldots,\widehat{r}_{NM-1}]^{T}$ of the
first row of the PN covariance matrix $\mathbf{R}(\tau)$ in (23) by using the
PN estimate $\widehat{\bm{\xi}}$:

$\displaystyle\widehat{r}_{i}=\frac{1}{NM-i}\sum_{k=0}^{NM-i-1}\widehat{\xi}_{k}\,\widehat{\xi}_{k+i}\leavevmode\nobreak\
,$ (59)

where
$\widehat{\bm{\xi}}=[\widehat{\xi}_{0},\ldots,\widehat{\xi}_{NM-1}]^{T}$.

  2. 2.

Construct the set of possible delay values by unfolding the ambiguous delay
estimate ${\widehat{\tau}}$ up to some maximum ambiguity index $K$:

$\displaystyle\mathcal{T}=\\{\tau\leavevmode\nobreak\
\lvert\leavevmode\nobreak\
\tau={\widehat{\tau}}+kT,\,k=0,\ldots,K\\}\leavevmode\nobreak\ .$ (60)

  3. 3.

Find the unambiguous (true) estimate of delay by solving the parametric
Toeplitz-block Toeplitz PN matrix reconstruction problem:

$\displaystyle{{\widehat{\tau}}^{\rm{true}}}=\arg\min_{\tau\in\mathcal{T}}\leavevmode\nobreak\
\big{\lVert}[\mathbf{R}(\tau)]_{0,:}-\widehat{\mathbf{r}}\big{\rVert}^{2}\leavevmode\nobreak\
.$ (61)

##  VI Simulation Results

In this section, we assess the performance of the proposed delay-Doppler
estimation algorithm for PN-impaired OFDM radar sensing using the mmWave
setting in Table I. In the observation model (11), the complex data symbols
$\mathbf{X}$ are chosen randomly from the QPSK constellation. To evaluate the
RMSE performances, we generate a total of $2500$ Monte Carlo realizations
consisting of every combination of $50$ independent realizations of PN vector
$\bm{\xi}$ and AWGN vector $\mathbf{z}$ in (26). Unless otherwise stated, we
consider a target with $R=30\,\rm{m}$ and $v=20\,\rm{m/s}$, and an oscillator
with $f_{\rm{3dB}}=200\,\rm{kHz}$ and $f_{\rm{loop}}=1\,\rm{MHz}$ (in the case
of PLL) [48]. Moreover, SNR is defined as
${\rm{SNR}}=\lvert\alpha\rvert^{2}/(2\sigma^{2})$ according to the model in
(11). To provide comparative performance analysis, we consider the following
benchmark processing schemes:

  * •

MAP-ISAA: The proposed hybrid ML/MAP estimation algorithm based on ISAA,
described in Algorithm 1.

  * •

2-D FFT: The standard 2-D FFT method used in OFDM radar processing [61, 60,
12, 40, 39], which corresponds to the optimal estimator in the ML sense in the
absence of PN121212The ML estimator in the absence of PN can be obtained as a
special case of the hybrid ML/MAP estimator in (38) for known PN matrix with
all elements being equal to $1$, i.e., $\mathbf{W}={{\bm{1}}}_{N\times M}$, as
mentioned in Sec. II-D. This special case has already been derived in
(52)–(54) for a given PN estimate $\widehat{\mathbf{W}}$. Inserting
$\widehat{\mathbf{W}}={{\bm{1}}}_{N\times M}$ and
$\widehat{\bm{\xi}}={{\bm{0}}}_{NM}$ into (54), one can readily derive the ML
estimator in (62), which can be implemented via 2-D FFT as $\mathbf{b}(\tau)$
in (9) and $\mathbf{c}(\nu)$ in (10) are DFT matrix columns for a uniformly
sampled delay-Doppler grid.:

$\displaystyle({\widehat{\tau}},{\widehat{\nu}})=\arg{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\max_{\tau,\nu}}\leavevmode\nobreak\
\Big{\lvert}\mathbf{b}^{H}(\tau)\left(\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\mathbf{Y}\right)\mathbf{c}(\nu)\Big{\rvert}^{2}\leavevmode\nobreak\
.$ (62)

  * •

2-D FFT (PN-free): The 2-D FFT method applied on the PN-free version of the
radar observations, $\mathbf{Y}^{\rm{PN-free}}$, in (12):

$\displaystyle({\widehat{\tau}},{\widehat{\nu}})=\arg\max_{\tau,\nu}\leavevmode\nobreak\
\Big{\lvert}\mathbf{b}^{H}(\tau)\left(\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\mathbf{Y}^{\rm{PN-
free}}\right)\mathbf{c}(\nu)\Big{\rvert}^{2}\leavevmode\nobreak\ .$ (63)

This benchmark will provide insights into PN-induced performance losses and
gains in delay-Doppler estimation.

Besides the above schemes, we also plot the hybrid CRBs [52, 53] to quantify
the theoretical performance bounds. The hybrid CRB in the presence and absence
of PN will be denoted as “CRB” and “CRB (PN-free)”, respectively, which
theoretically lower-bound the RMSE of “MAP-ISAA” and “2-D FFT (PN-free)”
algorithms.

Table I: OFDM Simulation Parameters Parameter | Value
---|---
Carrier Frequency, $f_{c}$ | $28\,\rm{GHz}$
Total Bandwidth, $B$ | $50\,\rm{MHz}$
Number of Subcarriers, $N$ | $256$
Number of Symbols, $M$ | $10$
Subcarrier Spacing, $\Delta f$ | $195.31\,\rm{kHz}$
Symbol Duration, $T$ | $5.12\,\mu\rm{s}$
Cyclic Prefix Duration, $T_{\rm{cp}}$ | $1.28\,\mu\rm{s}$
Total Symbol Duration, $T_{\rm{sym}}$ | $6.40\,\rm{\mu s}$

In what follows, we first evaluate the RMSE performances of the considered
processing schemes under various operating conditions with regard to SNR,
oscillator quality and target range. Then, we illustrate the convergence
behavior of the proposed algorithm in Algorithm 1. Finally, we demonstrate the
PN exploitation capability of the proposed approach in Algorithm 2.

###  VI-A Performance with respect to SNR

We first assess the performance of the proposed MAP-ISAA algorithm, along with
the FFT-based benchmark methods, with respect to SNR. Fig. 3 and Fig. 4 show
(for FRO and PLL architectures, respectively) the range and velocity RMSEs of
the considered schemes131313Due to discretization of search space and its
confinement to a finite interval for practical implementation of the
estimators, the RMSE might slightly fall below the CRB in certain rare
scenarios.. In terms of ranging performance in Fig. 3 and Fig. 4, the proposed
algorithm achieves the CRB and significantly outperforms the 2-D FFT benchmark
for both FRO and PLL architectures, especially at medium and high SNRs, where
an order-of-magnitude improvement in ranging accuracy can be observed. It is
seen that the performance of the standard FFT method saturates above a certain
SNR level, while the proposed approach can effectively utilize the prior
information on PN to compensate for its impact on the observations and avoid
such saturation behavior by attaining its theoretical lower bound (which
decreases with increasing SNR). Moreover, in compliance with the theoretical
bounds, MAP-ISAA exhibits ranging performance very close to that achieved in
the absence of PN, which evidences its remarkable PN compensation capability.

In contrast to their ranging performances, FRO and PLL display different
trends in velocity estimation, depicted in Fig. 3 and Fig. 4. For PLL, the
MAP-ISAA algorithm can provide noticeable improvements in velocity RMSE over
the FFT benchmark, leading to gains on the order of several $10\,\rm{cm/s}$,
while in case of FRO performance gains seem negligible with respect to the FFT
method. The same observation is also valid when comparing the CRB to the RMSE
of the FFT method since MAP-ISAA can get very close to the CRB for both types
of oscillators. This difference between FRO and PLL in velocity estimation can
be attributed to the following two facts: (i) velocity information in (11) is
extracted from slow-time (symbol-to-symbol) phase rotations represented by
$\mathbf{c}(\nu)$ in (10), and (ii) PN samples in different symbols are
uncorrelated for FRO, while they can be correlated for PLL, as shown in (23),
Lemma 2 and Fig. 2. Therefore, correlation of PN across symbols can be
exploited in PLL for accurate PN estimation and compensation in slow-time,
leading to better velocity estimates compared to FRO. In this respect, similar
ranging performances for FRO and PLL can be explained by pointing out high
correlation of PN in fast-time for both oscillator types, represented by the
block-diagonals $\mathbf{R}_{0}(\tau)$ in (23) and (25).

Comparing the asymptotic trends of the MAP-ISAA algorithm, as well as the
corresponding CRB, between range and velocity estimation in Fig. 3 and Fig. 4,
we observe plateau in velocity estimation performance as opposed to
monotonically decreasing errors for range estimation. In connection with this,
PN causes only slight degradation of ranging performance, whereas velocity
estimation can be severely degraded by PN. This is due to the fact that range
information is gathered from frequency-domain (or, equivalently, fast-time)
phase shifts across $\mathbf{b}(\tau)$ in (9), while velocity information
comes from slow-time phase shifts in (10). Since PN enjoys much higher
correlation in fast-time than in slow-time, such correlation can be utilized
to cancel out its effect in fast-time and accordingly provide accurate range
estimates.

Figure 3: 3 Range and 3 velocity RMSE with respect to SNR for FRO with
$f_{\rm{3dB}}=200\,\rm{kHz}$.

Figure 4: 4 Range and 4 velocity RMSE with respect to SNR for PLL with
$f_{\rm{3dB}}=200\,\rm{kHz}$ and $f_{\rm{loop}}=1\,\rm{MHz}$.

###  VI-B Performance with respect to Oscillator Quality

We now investigate the performance of the considered schemes with respect to
oscillator quality for both FRO and PLL implementations. To this end, in Fig.
5 we report accuracy in range, velocity and PN estimation141414The PN RMSE is
defined as
$\big{[}\sum_{c=1}^{C}\big{\lVert}\widehat{\bm{\xi}}_{c}-\bm{\xi}_{c}\big{\rVert}^{2}/(NMC)\big{]}^{1/2}$,
where $\widehat{\bm{\xi}}_{c}\in\mathbb{R}^{NM\times 1}$ and
$\bm{\xi}_{c}\in\mathbb{R}^{NM\times 1}$ denote, respectively, the estimated
and true PN at the ${c^{\rm{th}}}$ Monte Carlo realization and $C=2500$. For
the 2-D FFT method, we set $\widehat{\bm{\xi}}_{c}={{\bm{0}}}$. against the
$3\,\rm{dB}$ bandwidth $f_{\rm{3dB}}$ of FRO at ${\rm{SNR}}=20\,\rm{dB}$. It
can be observed that the ranging performance of the proposed MAP-ISAA
algorithm is highly robust in the face of worsening oscillator quality (i.e.,
increasing $f_{\rm{3dB}}$), maintaining almost the same RMSE level and
attaining the corresponding bounds over a wide range of $f_{\rm{3dB}}$ values.
On the other hand, the FFT benchmark suffers from a considerable performance
degradation as $f_{\rm{3dB}}$ increases due to increasing PN variance, shown
in Fig. 5. Additionally, range accuracy obtained by MAP-ISAA is very close to
that achievable in the absence of PN for all $f_{\rm{3dB}}$ values, which
demonstrates the effectiveness of the proposed PN estimation/compensation
approach in Algorithm 1. Regarding velocity RMSE, similar trends to those in
Sec. VI-A can be observed for both the CRB and the RMSE of MAP-ISAA, due to
lack of PN correlation in slow-time (i.e., the block-diagonal structure of PN
covariance matrix for FRO, as specified in Lemma 2).

Figure 5: 5 Range, 5 velocity and 5 PN RMSE with respect to $3\,\rm{dB}$
bandwidth $f_{\rm{3dB}}$ of FRO at ${\rm{SNR}}=20\,\rm{dB}$.

Next, we consider the RMSE performances for PLL synthesizers in Fig. 6, which
shows accuracies against the PLL loop bandwidth $f_{\rm{loop}}$. In terms of
both range and velocity RMSEs, the PN-ignorant FFT method performs worse with
decreasing $f_{\rm{loop}}$ since PN variance increases as $f_{\rm{loop}}$
decreases, as seen from (19) and Fig. 6. On the other hand, the proposed
algorithm attains the same ranging accuracy for the entire interval of
$f_{\rm{loop}}$ values from $10\,\rm{kHz}$ to $10\,\rm{MHz}$, proving its
robustness against the quality of PLL synthesizer. Hence, the performance gain
of MAP-ISAA in range estimation with respect to the FFT benchmark becomes more
pronounced with decreasing PLL quality (i.e., decreasing $f_{\rm{loop}}$).
Moreover, the ranging performance of the proposed approach under the impact of
PN is near that of the FFT benchmark that uses PN-free radar observations over
a broad range of $f_{\rm{loop}}$ values, which indicates almost perfect
cancellation of the effect of PN. Furthermore, the velocity RMSE curves in
Fig. 6 reveal that MAP-ISAA performs very close to CRB and can provide gains
on the order of $10\,\rm{cm/s}$ over the FFT benchmark. Finally, from Fig. 5
and Fig. 6, we note that MAP-ISAA can achieve the CRBs on range, velocity and
PN estimation for both FROs and PLLs under various levels of oscillator
quality.

Figure 6: 6 Range, 6 velocity and 6 PN RMSE with respect to loop bandwidth
$f_{\rm{loop}}$ of PLL with $f_{\rm{3dB}}=200\,\rm{kHz}$ at
${\rm{SNR}}=20\,\rm{dB}$.

###  VI-C Performance with respect to Target Range

To further highlight the benefits of the proposed algorithm, we explore the
impact of target range on the performance of the processing schemes under
consideration, recalling the delay-dependent statistics of PN derived in Sec.
III. In Fig. 7 and Fig. 8, we show the range-velocity RMSEs versus target
range for FRO and PLL, respectively. As can be noticed, the ranging
performance of the proposed approach is robust against increasing target range
for both FRO and PLL implementations, which suggest that Algorithm 1 can
successfully utilize the knowledge of delay-dependent PN covariance
$\mathbf{R}(\tau)$ to jointly estimate the coupled delay and PN parameters.
However, the FFT method experiences substantial loss in ranging accuracy as
target moves further away from radar, in compliance with monotonically
increasing variance of PN as a function of delay in (18) and (19). Moreover,
comparing MAP-ISAA and CRB against the PN-free benchmark, PN-induced
performance degradation is only marginal, which agrees with the observations
in Sec. VI-A and Sec. VI-B. Furthermore, looking at the velocity RMSEs in Fig.
7 and Fig. 8, we observe that improvements on the order of $\rm{m/s}$ can be
provided by the proposed algorithm over the FFT benchmark for PLL, whereas no
noticeable gains occur in case of FRO. This further corroborates the above-
mentioned insights on the impact of different PN correlation characteristics
of FRO and PLL (specified in (23) and (25)) onto velocity accuracy.

Figure 7: 7 Range, and 7 velocity RMSE with respect to target range for FRO
with $f_{\rm{3dB}}=150\,\rm{kHz}$ at ${\rm{SNR}}=10\,\rm{dB}$.

Figure 8: 8 Range, and 8 velocity RMSE with respect to target range for PLL
with $f_{\rm{3dB}}=150\,\rm{kHz}$ and $f_{\rm{loop}}=100\,\rm{kHz}$ at
${\rm{SNR}}=10\,\rm{dB}$.

###  VI-D Convergence Behavior of the Proposed Algorithm

To illustrate the convergence behavior of the proposed MAP-ISAA algorithm in
Algorithm 1 through iterated approximations, Fig. 9 demonstrates the evolution
of range, velocity and PN RMSEs over consecutive iterations for PLL. It is
observed that starting from the output of the FFT benchmark, Algorithm 1
monotonically converges to the corresponding CRBs on range, velocity and PN
estimation within few iterations, which proves the effectiveness of the
proposed iterated SAA approach. To explore the effect of iterated
approximations on PN tracking performance, in Fig. 10 we depict instances from
the PN process along with the PN estimates at the first and last iteration of
Algorithm 1 for FRO and PLL architectures151515Phase wrapping in PN
trajectories (i.e., phases outside $[-\pi,\ \pi]$) occurs very rarely and does
not affect the resulting RMSE performances.. Comparing the results of the
first and last iterations, we see that by applying iterated approximations in
Algorithm 1, PN tracking accuracy significantly improves around regions with
high fluctuations, which verifies the superiority of the proposed iterative
refinement approach.

Figure 9: Evolution of range, velocity and PN RMSEs through successive
iterations of the proposed algorithm at ${\rm{SNR}}=10\,\rm{dB}$ for PLL with
$f_{\rm{3dB}}=200\,\rm{kHz}$ and $f_{\rm{loop}}=1\,\rm{MHz}$. Initialized at
the FFT output, the proposed approach can quickly converge to the
corresponding CRBs in few iterations, recovering the performance loss incurred
by the PN.

Figure 10: True PN realization and PN estimates at the first and last
iteration of the proposed algorithm at ${\rm{SNR}}=10\,\rm{dB}$ for 10 FRO
with $f_{\rm{3dB}}=150\,\rm{kHz}$ considering a target with range
$R=100\,\rm{m}$, and 10 PLL with $f_{\rm{3dB}}=150\,\rm{kHz}$ and
$f_{\rm{loop}}=100\,\rm{kHz}$ considering a target with range $R=130\,\rm{m}$.
The proposed iterated SAA approach provides substantial improvements in PN
trajectory tracking performance especially when large fluctuations occur.

###  VI-E PN Exploitation Capability of the Proposed Algorithm

In this part, we demonstrate the PN exploitation capability of the proposed
method in Algorithm 2 by considering a target located at $R=1000\,\rm{m}$,
which is beyond the maximum unambiguous range of $768\,\rm{m}$ according to
Table I. In Fig. 11, we plot the range RMSE vs. SNR for PLL with
$f_{\rm{3dB}}=20\,\rm{kHz}$ and $f_{\rm{loop}}=1\,\rm{MHz}$ by searching for
the first two ambiguity intervals in (60), i.e., $K=1$. In addition, Fig. 11
shows the PN covariances at ${\rm{SNR}}=20\,\rm{dB}$, obtained via the
estimate $\widehat{\mathbf{r}}$ in (59) and via the analytical evaluation
$[\mathbf{R}(\tau)]_{0,:}$ in (61) at the ambiguous and true ranges. It can be
seen that the proposed PN exploitation approach can successfully resolve the
range ambiguity starting from ${\rm{SNR}}=15\,\rm{dB}$ (when the PN estimate
becomes sufficiently accurate161616Additional simulations with different
$f_{\rm{3dB}}$ values for PLL offer an intriguing insight into the relation
between $f_{\rm{3dB}}$ and the degree of accuracy improvement via PN
exploitation: higher $f_{\rm{3dB}}$ (which means lower oscillator quality)
provides improved resolvability of range ambiguity at low SNRs through more
pronounced notches in the PN covariance profile at the ambiguous and true
ranges, while leading to lower accuracy at high SNRs due to performance
saturation (i.e., when the ambiguity is already resolved).) and attain the
corresponding CRB, whereas neither of the FFT benchmarks can correctly
identify the true target range due to intrinsic range ambiguity in
$\mathbf{b}(\tau)$, in compliance with the discussions in Sec. V. More
specifically, the FFT processing in (62), either using PN-impaired
observations $\mathbf{Y}$ in (11) or PN-free observations $\mathbf{Y}^{\rm{PN-
free}}$ in (12), can extract the range information only from
$\mathbf{b}(\tau)$, which causes ambiguity for $\tau>T$ (see (9)). On the
other hand, to estimate the true (unambiguous) range of the target, Algorithm
2 can effectively exploit the ranging information revealed by the PN
statistics $\mathbf{R}(\tau)$ in (20), which does not introduce any ambiguity
in estimating $\tau$ as noticed from (17)–(21). Surprisingly, the proposed
approach under the impact of PN can significantly outperform the FFT benchmark
that uses ideal, PN-free observations, indicating its excellent PN
exploitation performance.

Figure 11: PN exploitation to resolve range ambiguity for PLL with
$f_{\rm{3dB}}=20\,\rm{kHz}$ and $f_{\rm{loop}}=1\,\rm{MHz}$, where target is
located at $R=1000\,\rm{m}$. 11 Range RMSE with respect to SNR. 11 PN
covariance at ${\rm{SNR}}=20\,\rm{dB}$, corresponding to the estimate
$\widehat{\mathbf{r}}$ in (59) and the analytical expression
$[\mathbf{R}(\tau)]_{0,:}$ in (61) evaluated at $\tau_{1}={\widehat{\tau}}$
and $\tau_{2}={\widehat{\tau}}+T$, where ${\widehat{\tau}}$ is the ambiguous
delay estimate. The proposed approach in Algorithm 2 can correctly identify
the true target range by exploiting PN statistics ($\widehat{\mathbf{r}}$
matches with $[\mathbf{R}(\tau_{2})]_{0,:}$), while the FFT benchmark fails to
resolve range ambiguity irrespective of whether the PN is present or absent.

##  VII Concluding Remarks

In this work, we have investigated the problem of monostatic radar sensing in
OFDM JRC systems under the effect of PN. Starting from an explicit derivation
of PN statistics in the OFDM radar receiver, we have proposed a novel
algorithm for joint estimation of delay, Doppler and PN by devising an
iterated small angle approximation approach that solves the hybrid ML/MAP
optimization problem through alternating updates of delay-Doppler and PN
estimates. In addition, we have developed a PN exploitation algorithm that
uses the statistics of the PN estimates to resolve range ambiguity. To assess
the performance of the proposed algorithms, comprehensive simulations have
been carried out over a broad range of operating conditions, leading to the
following key findings:

  * •

Range vs. Velocity Accuracy under PN: The proposed MAP-ISAA algorithm
significantly outperforms the FFT benchmark in ranging accuracy, while
providing slight improvements in velocity accuracy. This results from strong
(weak) correlation of PN in fast-time (slow-time).

  * •

Impact of Oscillator Type - FRO vs. PLL: FROs and PLLs lead to very similar
performance in range estimation, whereas PLLs can provide much higher velocity
accuracy than FROs. This is due to the absence of PN correlation in slow-time
for FROs.

  * •

Scenarios with Significant Impact of PN on Sensing: The effect of PN on
sensing performance becomes more significant at higher SNRs, for oscillators
with larger $f_{\rm{3dB}}$ and smaller $f_{\rm{loop}}$ values, and farther
targets. In such scenarios, PN should be considered in algorithm design.

  * •

Iterated Small Angle Approximation: The proposed algorithm enjoys fast
convergence to the hybrid CRB and provides superior PN tracking performance
via iterated approximations.

  * •

PN Exploitation: Delay-dependency of PN statistics can be effectively
exploited to resolve range ambiguity above a certain SNR level, thereby
converting the detrimental effect of PN into something beneficial for radar
sensing.

As future research, we plan to investigate extensions of the proposed
framework to MIMO architectures and the multi-target case.

##  Appendix A Correlation Function of Differential Phase Noise in (17)

Let the correlation function of $\phi(t)$ in (14) be defined as

$\displaystyle R_{\phi\phi}(t_{1},t_{2})$
$\displaystyle\triangleq\operatorname*{\mathbb{E}}\left[\phi(t_{1})\phi(t_{2})\right]\leavevmode\nobreak\
.$ (64)

Then, it follows from (15) that

$\displaystyle\operatorname*{\mathbb{E}}\left[\xi(t_{1},\tau)\xi(t_{2},\tau)\right]$
$\displaystyle=R_{\phi\phi}(t_{1},t_{2})-R_{\phi\phi}(t_{1},t_{2}-\tau)$ (65)
$\displaystyle-
R_{\phi\phi}(t_{1}-\tau,t_{2})+R_{\phi\phi}(t_{1}-\tau,t_{2}-\tau)\leavevmode\nobreak\
.$

From (15) and (16), we have

$\displaystyle{\sigma_{\xi}^{2}(\tau)}=\operatorname*{\mathbb{E}}\left[\xi(t,\tau)^{2}\right]=\sigma^{2}_{\phi}(t)+\sigma^{2}_{\phi}(t-\tau)-2R_{\phi\phi}(t,t-\tau)\leavevmode\nobreak\
,$ (66)

which yields

$\displaystyle
R_{\phi\phi}(t_{1},t_{2})=\frac{\sigma^{2}_{\phi}(t_{1})+\sigma^{2}_{\phi}(t_{2})-\sigma^{2}_{\xi}(t_{1}-t_{2})}{2}\leavevmode\nobreak\
.$ (67)

Finally, inserting (67) into (65) yields

$\displaystyle\operatorname*{\mathbb{E}}\left[\xi(t_{1},\tau)\xi(t_{2},\tau)\right]$
(68)
$\displaystyle=\frac{\sigma^{2}_{\xi}(t_{1}-t_{2}+\tau)+\sigma^{2}_{\xi}(t_{1}-t_{2}-\tau)}{2}-\sigma^{2}_{\xi}(t_{1}-t_{2})\leavevmode\nobreak\
.$

Defining $\Delta t\triangleq t_{1}-t_{2}$ in (68) and considering
$\sigma^{2}_{\xi}(\tau)=\sigma^{2}_{\xi}(-\tau)$, we obtain the result in
(17).

##  Appendix B Block-Diagonal Structure of $\mathbf{R}(\tau)$ for FROs

Using (18) in (17), we obtain the correlation function for FROs as

$\displaystyle R_{\xi\xi}(\Delta t,\tau)$ $\displaystyle=4\pi
f_{\rm{3dB}}\max\left(\tau-\lvert\Delta t\rvert,0\right)\leavevmode\nobreak\
.$ (69)

Since $\lvert\Delta t^{(m)}_{n_{1}n_{2}}\rvert\geq T_{\rm{cp}}$ for $m>0$ in
(24), we have $R_{\xi\xi}(\Delta t^{(m)}_{n_{1}n_{2}},\tau)=0$ for $m>0$ and
$\tau\leq T_{\rm{cp}}$ from (69), leading to $\mathbf{R}_{m}(\tau)={{\bm{0}}}$
for $m>0$, which yields (25).

##  Appendix C Obtaining (48) From (47)

In this part, we provide details on how to reach (48) from (47). To this end,
we present the following auxiliary lemmas.

######  Lemma 3.

For any $\mathbf{x}\in\mathbb{R}^{N\times 1}$ and any Hermitian matrix
$\mathbf{A}\in\mathbb{C}^{N\times N}$, the following equality holds:

$\displaystyle\mathbf{x}^{T}\mathbf{A}\mathbf{x}=\mathbf{x}^{T}\Re\left\\{\mathbf{A}\right\\}\mathbf{x}\leavevmode\nobreak\
.$ (70)

######  Proof.

Since $\mathbf{x}^{T}\mathbf{A}\mathbf{x}$ is scalar, we can write

$\displaystyle\mathbf{x}^{T}\mathbf{A}\mathbf{x}=(\mathbf{x}^{T}\mathbf{A}\mathbf{x})^{T}=\mathbf{x}^{T}\mathbf{A}^{T}\mathbf{x}\stackrel{{\scriptstyle\textnormal{(a)}}}{{\mathstrut{=}}}\mathbf{x}^{T}\mathbf{A}^{\ast}\mathbf{x}\leavevmode\nobreak\
,$ (71)

where (71) is due to $\mathbf{A}$ being Hermitian, i.e.,
$\mathbf{A}^{H}=\mathbf{A}$. Re-writing the leftmost and rightmost terms in
(71), we obtain

$\displaystyle\mathbf{x}^{T}\mathbf{A}\mathbf{x}=\mathbf{x}^{T}\left(\Re\left\\{\mathbf{A}\right\\}+j\Im\left\\{\mathbf{A}\right\\}\right)\mathbf{x}$
(72)
$\displaystyle=\mathbf{x}^{T}\mathbf{A}^{\ast}\mathbf{x}=\mathbf{x}^{T}\left(\Re\left\\{\mathbf{A}\right\\}-j\Im\left\\{\mathbf{A}\right\\}\right)\mathbf{x}\leavevmode\nobreak\
,$

which leads to

$\displaystyle\mathbf{x}^{T}\Im\left\\{\mathbf{A}\right\\}\mathbf{x}=0\leavevmode\nobreak\
.$ (73)

Substituting (73) into (72) yields the result in (74). $\blacksquare$

######  Lemma 4.

For any $\mathbf{x}\in\mathbb{R}^{N\times 1}$,
$\mathbf{y}\in\mathbb{R}^{N\times 1}$ and any Hermitian matrix
$\mathbf{A}\in\mathbb{C}^{N\times N}$, the following equality holds:

$\displaystyle
j(\mathbf{x}^{T}\mathbf{A}\mathbf{y}-\mathbf{y}^{T}\mathbf{A}\mathbf{x})=2\mathbf{y}^{T}\Im\left\\{\mathbf{A}\right\\}\mathbf{x}\leavevmode\nobreak\
.$ (74)

######  Proof.

Since $\mathbf{x}^{T}\mathbf{A}\mathbf{y}$ is scalar and $\mathbf{A}$ is
Hermitian, we can write

$\displaystyle
j(\mathbf{x}^{T}\mathbf{A}\mathbf{y}-\mathbf{y}^{T}\mathbf{A}\mathbf{x})=j(\mathbf{y}^{T}\mathbf{A}^{T}\mathbf{x}-\mathbf{y}^{T}\mathbf{A}\mathbf{x})$
$\displaystyle=j(\mathbf{y}^{T}\mathbf{A}^{\ast}\mathbf{x}-\mathbf{y}^{T}\mathbf{A}\mathbf{x})$
$\displaystyle=j\mathbf{y}^{T}(\mathbf{A}^{\ast}-\mathbf{A})\mathbf{x}$
$\displaystyle=-j\mathbf{y}^{T}2j\Im\left\\{\mathbf{A}\right\\}\mathbf{x}$
$\displaystyle=2\mathbf{y}^{T}\Im\left\\{\mathbf{A}\right\\}\mathbf{x}\leavevmode\nobreak\
,$

which completes the proof. $\blacksquare$

Opening up the first term in (47), we have

$\displaystyle({{\bm{1}}}-j\bm{\xi}_{\Delta})^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)({{\bm{1}}}+j\bm{\xi}_{\Delta})$
(75)
$\displaystyle={{\bm{1}}}^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu){{\bm{1}}}+\bm{\xi}_{\Delta}^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\bm{\xi}_{\Delta}$
(76) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+j\left({{\bm{1}}}^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\bm{\xi}_{\Delta}-\bm{\xi}_{\Delta}^{T}\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu){{\bm{1}}}\right)$
$\displaystyle={{\bm{1}}}^{T}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+\bm{\xi}_{\Delta}^{T}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\bm{\xi}_{\Delta}$
(77) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+2\bm{\xi}_{\Delta}^{T}\Im\left\\{\mathbf{\Gamma}(\tau,\nu\right\\}{{\bm{1}}}\leavevmode\nobreak\
,$

where the transition from (76) to (77) is via Lemma 3 and Lemma 4. Plugging
(75)–(77) into (47) yields (48).

##  Appendix D Delay-Doppler Estimation in (52)

Using (27), (35) and (39) in (52), we obtain

$\displaystyle\mathcal{L}(\tau,\nu,\widehat{\bm{\xi}})$
$\displaystyle=\frac{1}{\sigma^{2}}\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)^{H}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)+\widehat{\bm{\xi}}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\log\det\mathbf{R}(\tau)$
$\displaystyle=\frac{1}{\sigma^{2}}\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)^{H}\left({\bm{\mathrm{I}}}-\frac{\mathbf{q}(\tau,\nu)\mathbf{q}^{H}(\tau,\nu)}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\right)\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\widehat{\bm{\xi}}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}+\log\det\mathbf{R}(\tau)$
$\displaystyle=\frac{\left\lVert\mathbf{y}\right\rVert^{2}}{\sigma^{2}}-\frac{\Big{\lvert}\mathbf{q}^{H}(\tau,\nu)\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)\Big{\rvert}^{2}}{\sigma^{2}\left\lVert\mathbf{X}\right\rVert_{F}^{2}}+\widehat{\bm{\xi}}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\log\det\mathbf{R}(\tau)\leavevmode\nobreak\ .$ (78)

Based on (78), the problem (52) is equivalent to

$\displaystyle({\widehat{\tau}},{\widehat{\nu}})$
$\displaystyle=\arg\max_{\tau,\nu}\leavevmode\nobreak\
\frac{\Big{\lvert}\mathbf{q}^{H}(\tau,\nu)\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)\Big{\rvert}^{2}}{\sigma^{2}\left\lVert\mathbf{X}\right\rVert_{F}^{2}}-\widehat{\bm{\xi}}^{T}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-\log\det\mathbf{R}(\tau)\leavevmode\nobreak\ .$ (79)

Using (28) in the first term of (79) and defining
$\widehat{\mathbf{W}}\triangleq{\rm{reshape}}_{N,M}\left(e^{-j\widehat{\bm{\xi}}}\right)$,
we have

$\displaystyle\Big{\lvert}\mathbf{q}^{H}(\tau,\nu)\left(e^{j\widehat{\bm{\xi}}}\odot\mathbf{y}\right)\Big{\rvert}^{2}$
(80)
$\displaystyle=\Big{\lvert}{{{\rm{tr}}\bigg{(}\left[\mathbf{F}_{N}^{H}\Big{(}\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\Big{)}\right]^{H}\left[\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\right]\bigg{)}}}\Big{\rvert}^{2}$
$\displaystyle=\Big{\lvert}{{{\rm{tr}}\bigg{(}\Big{(}\mathbf{X}\odot\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\Big{)}^{H}\mathbf{F}_{N}\left[\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\right]\bigg{)}}}\Big{\rvert}^{2}$
$\displaystyle=\Big{\lvert}{{{\rm{tr}}\bigg{(}\Big{(}\mathbf{b}(\tau)\mathbf{c}^{H}(\nu)\Big{)}^{H}\Big{(}\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\left[\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\right]\Big{)}\bigg{)}}}\Big{\rvert}^{2}$
$\displaystyle=\Big{\lvert}{{{\rm{tr}}\bigg{(}\mathbf{c}(\nu)\mathbf{b}^{H}(\tau)\Big{(}\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\left[\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\right]\Big{)}\bigg{)}}}\Big{\rvert}^{2}$
$\displaystyle=\Big{\lvert}{{{\rm{tr}}\bigg{(}\mathbf{b}^{H}(\tau)\Big{(}\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\left[\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\right]\Big{)}\mathbf{c}(\nu)\bigg{)}}}\Big{\rvert}^{2}$
$\displaystyle=\Big{\lvert}\mathbf{b}^{H}(\tau)\left(\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\big{(}\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\big{)}\right)\mathbf{c}(\nu)\Big{\rvert}^{2}\leavevmode\nobreak\
,$

which, after being inserted into (79), results in (53).

##  Appendix E Convergence Analysis of Algorithm 1

In this part, we provide the convergence analysis of the proposed ISAA
algorithm in Algorithm 1 to identify the conditions under which the algorithm
converges to a stationary point of the hybrid ML/MAP optimization problem in
(38).

###  E-A Formulation of Subproblems to Solve (38)

In Algorithm 1, we decompose the original problem (38) into two subproblems,
each optimizing over either delay-Doppler or PN while keeping the other
variable fixed in an alternating fashion. Let
$({\widehat{\tau}}^{(i)},{\widehat{\nu}}^{(i)},\widehat{\bm{\xi}}^{(i)})$
denote the delay, Doppler and PN estimates at the beginning of the
${i^{\rm{th}}}$ iteration of the alternating optimization procedure.

####  E-A1 Subproblem for PN Estimation

The subproblem of (38) for PN estimation at the ${i^{\rm{th}}}$ iteration can
be cast as

$\displaystyle\widehat{\bm{\xi}}^{(i+1)}=\arg\min_{\bm{\xi}}\leavevmode\nobreak\
\mathcal{L}({\widehat{\tau}}^{(i)},{\widehat{\nu}}^{(i)},\bm{\xi})\leavevmode\nobreak\
,$ (81)

where $\mathcal{L}(\tau,\nu,\bm{\xi})$ is the hybrid ML/MAP objective function
in (39). In the proposed ISAA approach, we tackle (81) by solving the
equivalent problem

$\displaystyle\widehat{\bm{\xi}}_{\Delta}=\arg\min_{\bm{\xi}_{\Delta}}\leavevmode\nobreak\
\mathcal{L}({\widehat{\tau}}^{(i)},{\widehat{\nu}}^{(i)},\bm{\xi}_{\Delta}+\widehat{\bm{\xi}}^{(i)})\leavevmode\nobreak\
,$ (82)

which estimates the residual PN
$\bm{\xi}_{\Delta}=\bm{\xi}-\widehat{\bm{\xi}}^{(i)}$ in (40) instead of the
actual PN $\bm{\xi}$ given the PN estimate $\widehat{\bm{\xi}}^{(i)}$ from the
previous iteration. In the PN update step in Line 6 of Algorithm 1, the
problem (81), or, equivalently (82), is solved in closed-form by applying the
SAA for $\bm{\xi}_{\Delta}$:

$\displaystyle
e^{-j\bm{\xi}_{\Delta}}\approx{{\bm{1}}}-j\bm{\xi}_{\Delta}\leavevmode\nobreak\
.$ (83)

####  E-A2 Subproblem for Delay-Doppler Estimation

The subproblem of (38) for delay-Doppler estimation at the ${i^{\rm{th}}}$
iteration is given by

$\displaystyle({\widehat{\tau}}^{(i+1)},{\widehat{\nu}}^{(i+1)})=\arg\min_{\tau,\nu}\leavevmode\nobreak\
\mathcal{L}(\tau,\nu,\widehat{\bm{\xi}}^{(i+1)})\leavevmode\nobreak\ .$ (84)

The delay-Doppler update in Line 7 of Algorithm 1 finds an optimal solution of
the problem (84).

###  E-B Convergence Result

Based on the subproblem definitions in Sec. E-A, we can present the
convergence result for Algorithm 1.

######  Proposition 1.

Assume that the SAA in (83) holds (i.e., the residual PN is small) and there
exists a unique delay-Doppler pair that minimizes the objective in (84). Then,
Algorithm 1, which solves (81) and (84) in an alternating fashion, converges
to a stationary point of the hybrid ML/MAP optimization problem in (38).

######  Proof.

When the SAA in (83) holds, the closed-form solution given by (50) and (51)
provides the optimal solution of the PN estimation subproblem in (81), or,
equivalently (82), and is therefore a unique solution (being a closed-form
one). Then, relying on the assumption that the solution to (84) is unique, it
follows from [68, Proposition 1] that Algorithm 1 converges to a stationary
point of (38).

##  Appendix F Complexity Analysis of Algorithm 1

In this section, we analyze the per-iteration complexity of Algorithm 1. We
first focus on the update of PN estimates via (50) and (51), and then study
the complexity of updating delay-Doppler estimates via (53) and (54).

###  F-A Complexity of PN Estimation in (50)–(51)

For convenience, we repeat the residual PN estimate in (50) here:

$\displaystyle\widehat{\bm{\xi}}_{\Delta}(\tau,\nu)$
$\displaystyle=-\mathbf{R}(\tau)\Big{(}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)+\sigma^{2}{\bm{\mathrm{I}}}\Big{)}^{-1}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\times\left(\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+\sigma^{2}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}\right)\leavevmode\nobreak\
.$ (85)

In the following, the complexity of (85) is analyzed in four steps.

####  F-A1 Complexity of $\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$

For efficient computation of $\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$,
the conjugate gradient (CG) method can be employed, similar to [29]. Note that
computing

$\displaystyle\mathbf{x}=\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$ (86)

is equivalent to solving the following linear system of equations for
$\mathbf{x}$ [29]:

$\displaystyle\mathbf{R}(\tau)\mathbf{x}=\widehat{\bm{\xi}}^{(i)}\leavevmode\nobreak\
.$ (87)

The CG method provides an iterative procedure to solve (87), where at each
iteration, the major complexity results from a matrix-vector multiplication
[29]

$\displaystyle\bm{\omega}=\mathbf{R}(\tau)\bm{\kappa}$ (88)

for some $\bm{\kappa}\in\mathbb{R}^{NM\times 1}$, where
$\mathbf{R}(\tau)\in\mathbb{R}^{NM\times NM}$ is the PN covariance matrix
given in (23).

Let us define
$\bm{\omega}=[\bm{\omega}_{0}^{T}\,\ldots\,\bm{\omega}_{M-1}^{T}]^{T}$ and
$\bm{\kappa}=[\bm{\kappa}_{0}^{T}\,\ldots\,\bm{\kappa}_{M-1}^{T}]^{T}$, where
$\bm{\omega}_{m}\in\mathbb{R}^{N\times 1}$ and
$\bm{\kappa}_{m}\in\mathbb{R}^{N\times 1}$ for each $m$. Then, using the
Toeplitz-block Toeplitz structure of $\mathbf{R}(\tau)$ in (23), we can re-
write (88) as

$\displaystyle\bm{\omega}_{m}=\sum_{\ell=0}^{m-1}\mathbf{R}_{m-\ell}^{T}(\tau)\bm{\kappa}_{\ell}+\sum_{\ell=m}^{M-1}\mathbf{R}_{\ell-m}(\tau)\bm{\kappa}_{\ell}$
(89)

for $m=0,\ldots,M-1$, where $\mathbf{R}_{m}(\tau)\in\mathbb{R}^{N\times N}$ is
a Toeplitz matrix. Using the circulant approximation of Toeplitz matrices for
large $N$ [69], $\mathbf{R}_{m}(\tau)$ can be approximated as

$\displaystyle\mathbf{R}_{m}(\tau)\approx\mathbf{R}^{{\rm{circ}}}_{m}(\tau)=\mathbf{F}_{N}{\bm{\Lambda}}_{m}(\tau)\mathbf{F}_{N}^{H}\leavevmode\nobreak\
,$ (90)

where $\mathbf{F}_{N}\in\mathbb{C}^{N\times N}$ is the unitary DFT matrix,
$\mathbf{R}^{{\rm{circ}}}_{m}(\tau)\in\mathbb{R}^{N\times N}$ is the best
circulant approximation to $\mathbf{R}_{m}(\tau)$ that minimizes
$\left\lVert\mathbf{R}_{m}(\tau)-\mathbf{R}^{{\rm{circ}}}_{m}(\tau)\right\rVert_{F}$
[69], and

$\displaystyle{\bm{\Lambda}}_{m}(\tau)={\rm{diag}}\left(\sqrt{N}\mathbf{F}_{N}^{H}\mathbf{c}_{m}(\tau)\right)\in\mathbb{C}^{N\times
N}\leavevmode\nobreak\ ,$ (91)

with $\mathbf{c}_{m}(\tau)$ denoting the first column of
$\mathbf{R}^{{\rm{circ}}}_{m}(\tau)$.

Based on the approximation in (90), each summand in the second term of (89)
can be expressed as

$\displaystyle\mathbf{R}_{\ell-m}(\tau)\bm{\kappa}_{\ell}\approx\mathbf{F}_{N}{\bm{\Lambda}}_{\ell-m}(\tau)\mathbf{F}_{N}^{H}\bm{\kappa}_{\ell}\leavevmode\nobreak\
,$ (92)

which can be computed efficiently using FFT and IFFT, leading to
$\mathcal{O}(N\log N)$ complexity. Similarly, for the summands in the first
term of (89), we must compute

$\displaystyle\mathbf{R}_{m-\ell}^{T}(\tau)\bm{\kappa}_{\ell}\approx\mathbf{F}_{N}^{H}{\bm{\Lambda}}_{m-\ell}(\tau)\mathbf{F}_{N}\bm{\kappa}_{\ell}\leavevmode\nobreak\
,$ (93)

which results in $\mathcal{O}(N\log N)$ complexity. In practice, as seen from
the correlation characteristics of PN in Fig. 2, the PN correlation vanishes
after a certain time lag, meaning that only a few blocks in
$\mathbf{R}(\tau)$, say $M_{0}$, are dominant in the computation of the right-
hand side of (89) (i.e.,
$\mathbf{R}_{0}(\tau),\ldots,\mathbf{R}_{M_{0}-1}(\tau)$), where $M_{0}=1$ for
FROs due to the block-diagonal structure derived in Lemma 2 and $1\leq
M_{0}\ll M$ for PLLs (typically, $M_{0}\leq 3$). Hence, the complexity of (89)
is given by

$\displaystyle\mathcal{O}(M_{0}N\log N)\leavevmode\nobreak\ ,$ (94)

which yields the following complexity for (88):

$\displaystyle\mathcal{O}(M_{0}MN\log N)\leavevmode\nobreak\ .$ (95)

Assuming that the CG method requires $I$ iterations to converge, the
complexity of evaluating $\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$ in
(86) can be expressed as

$\displaystyle\mathcal{O}(IM_{0}MN\log N)\leavevmode\nobreak\ .$ (96)

####  F-A2 Complexity of
$\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}$

Using the definition of
$\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)$
in (45), we can write

$\displaystyle\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)$
(97)
$\displaystyle=\left({\rm{diag}}\left(\mathbf{y}\right)^{H}\bm{\Pi}^{\perp}_{\mathbf{q}(\tau,\nu)}{\rm{diag}}\left(\mathbf{y}\right)\right)\odot\left(e^{-j\widehat{\bm{\xi}}^{(i)}}(e^{j\widehat{\bm{\xi}}^{(i)}})^{T}\right)$
$\displaystyle=\left[{\rm{diag}}\left(\mathbf{y}\right)^{H}{\rm{diag}}\left(\mathbf{y}\right)-\left(\mathbf{y}^{\ast}\mathbf{y}^{T}\odot\frac{\mathbf{q}(\tau,\nu)\mathbf{q}^{H}(\tau,\nu)}{\left\lVert\mathbf{q}(\tau,\nu)\right\rVert^{2}}\right)\right]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\odot\left(e^{-j\widehat{\bm{\xi}}^{(i)}}(e^{j\widehat{\bm{\xi}}^{(i)}})^{T}\right)$
$\displaystyle={\rm{diag}}\left(\lvert\mathbf{y}\rvert^{2}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{1}{\left\lVert\mathbf{q}(\tau,\nu)\right\rVert^{2}}\left(\mathbf{y}^{\ast}\odot\mathbf{q}(\tau,\nu)\odot
e^{-j\widehat{\bm{\xi}}^{(i)}}\right)\left(\mathbf{y}^{\ast}\odot\mathbf{q}(\tau,\nu)\odot
e^{-j\widehat{\bm{\xi}}^{(i)}}\right)^{H}$
$\displaystyle={\rm{diag}}\left(\lvert\mathbf{y}\rvert^{2}\right)-\frac{1}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)^{H}$

where
$\widehat{\mathbf{\Xi}}^{(i)}\triangleq{\rm{diag}}\left(e^{-j\widehat{\bm{\xi}}^{(i)}}\right)$.
Using (97), the real and imaginary parts of
$\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)$
are given by

$\displaystyle\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}={\rm{diag}}\left(\lvert\mathbf{y}\rvert^{2}\right)$
(98) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{1}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\Re\left\\{\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)^{H}\right\\}\leavevmode\nobreak\
,$
$\displaystyle\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}$
(99) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
=-\frac{1}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\Im\left\\{\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)^{H}\right\\}\leavevmode\nobreak\
.$

Using (99), we can compute
$\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}$
as

$\displaystyle\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}$
(100)
$\displaystyle=-\frac{1}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\Im\left\\{\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)^{H}{{\bm{1}}}\right\\}\leavevmode\nobreak\
,$

which requires

$\displaystyle\mathcal{O}(NM)$ (101)

operations.

####  F-A3 Complexity of
$\Big{(}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)+\sigma^{2}{\bm{\mathrm{I}}}\Big{)}^{-1}\bm{\varpi}$

Now that we have computed $\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}$ and
$\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}$
in (85), we will analyze the complexity of evaluating

$\displaystyle\Big{(}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)+\sigma^{2}{\bm{\mathrm{I}}}\Big{)}^{-1}\bm{\varpi}\leavevmode\nobreak\
,$ (102)

where

$\displaystyle\bm{\varpi}=\Im\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}{{\bm{1}}}+\sigma^{2}\mathbf{R}(\tau)^{-1}\widehat{\bm{\xi}}^{(i)}\leavevmode\nobreak\
.$ (103)

Similar to the above analysis, the CG method can be employed to evaluate
(102), where each iteration involves the computation of

$\displaystyle\Big{(}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)+\sigma^{2}{\bm{\mathrm{I}}}\Big{)}\bm{\kappa}$
(104)

for some $\bm{\kappa}\in\mathbb{R}^{NM\times 1}$. Opening up the terms in
(104), we have

$\displaystyle\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)\bm{\kappa}+\sigma^{2}\bm{\kappa}\leavevmode\nobreak\
.$ (105)

As previously derived in (95), evaluation of $\mathbf{R}(\tau)\bm{\kappa}$ in
(105) has the complexity of

$\displaystyle\mathcal{O}(M_{0}MN\log N)\leavevmode\nobreak\ .$ (106)

Next, defining $\bm{\omega}=\mathbf{R}(\tau)\bm{\kappa}$, we compute
$\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\bm{\omega}$
in (105) using (98) as

$\displaystyle\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\bm{\omega}$
(107)
$\displaystyle={\rm{diag}}\left(\lvert\mathbf{y}\rvert^{2}\right)\bm{\omega}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{1}{\left\lVert\mathbf{X}\right\rVert_{F}^{2}}\Re\left\\{\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)\left(\mathbf{y}^{\ast}\odot\widehat{\mathbf{\Xi}}^{(i)}\mathbf{q}(\tau,\nu)\right)^{H}\bm{\omega}\right\\}\leavevmode\nobreak\
,$

which leads to a complexity of $\mathcal{O}(NM)$. Combining this with (106),
the complexity of computing the expression in (104) becomes

$\displaystyle\mathcal{O}(M_{0}MN\log N)\leavevmode\nobreak\ .$ (108)

Finally, assuming $I$ iterations for the CG method to converge, (102) has a
computational complexity of

$\displaystyle\mathcal{O}(IM_{0}MN\log N)\leavevmode\nobreak\ .$ (109)

####  F-A4 Final Complexity Result for (85)

After evaluating (102), the final step to obtain (85) is to compute

$\displaystyle\mathbf{R}(\tau)\bm{\varsigma}\leavevmode\nobreak\ ,$ (110)

where

$\displaystyle\bm{\varsigma}=\Big{(}\Re\left\\{\mathbf{\Gamma}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(i)}}(\tau,\nu)\right\\}\mathbf{R}(\tau)+\sigma^{2}{\bm{\mathrm{I}}}\Big{)}^{-1}\bm{\varpi}\leavevmode\nobreak\
,$ (111)

with $\bm{\varpi}$ being defined in (103). From (95), evaluation of (110)
involves

$\displaystyle\mathcal{O}(M_{0}MN\log N)$ (112)

operations. Combining (96), (101), (109) and (112), the computational
complexity of (85) is obtained as

$\displaystyle\mathcal{O}(IM_{0}MN\log N)\leavevmode\nobreak\ .$ (113)

###  F-B Complexity of Delay-Doppler Estimation in (53)–(54)

To evaluate (54) in practice, it is enough to compute the first term as the
remaining terms are negligibly small and can be discarded. We notice that the
first term in (54) (up to a multiplicative constant)

$\displaystyle\Big{\lvert}\mathbf{b}^{H}(\tau)\left(\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\big{(}\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\big{)}\right)\mathbf{c}(\nu)\Big{\rvert}^{2}$
(114)

can be evaluated over a delay-Doppler grid via

$\displaystyle\Big{\lvert}\mathbf{F}_{N}^{H}\left(\mathbf{X}^{\ast}\odot\mathbf{F}_{N}\big{(}\widehat{\mathbf{W}}^{\ast}\odot\mathbf{Y}\big{)}\right)\mathbf{F}_{M}\Big{\rvert}^{2}\leavevmode\nobreak\
,$ (115)

utilizing the fact that $\mathbf{b}(\tau)$ in (9) and $\mathbf{c}(\nu)$ in
(10) correspond to DFT matrix columns for uniformly sampled delay-Doppler
values. The complexity of (115) is thus given by

$\displaystyle\mathcal{O}(NM)+\mathcal{O}(MN\log N)+\mathcal{O}(NM\log
M)\leavevmode\nobreak\ ,$ (116)

where $\mathcal{O}(NM)$ is due to element-wise multiplication of $N\times M$
matrices, $\mathcal{O}(MN\log N)$ results from $M$ times $N$-point FFT/IFFT,
and $\mathcal{O}(NM\log M)$ comes from $N$ times $M$-point FFT. Re-writing
(116), the complexity of (115) can be expressed as

$\displaystyle\mathcal{O}(NM\log(NM))\leavevmode\nobreak\ .$ (117)

###  F-C Complexity of Algorithm 1

Using the complexity of PN estimation in (113) and that of delay-Doppler
estimation in (117), the per-iteration complexity of Algorithm 1 can be
written as

$\displaystyle\mathcal{O}\big{(}MN\left(\log M+(IM_{0}+1)\log
N\right)\big{)}\leavevmode\nobreak\ .$ (118)

##  Appendix G Theoretical Ground for Formulation of PN Exploitation Problem
in (56)

In this section, we provide the theoretical motivation behind the parametric
covariance matrix reconstruction problem in (56), which is formulated to
resolve range ambiguity via PN exploitation. The goal of PN exploitation is to
estimate the unknown delay $\tau$ from the PN estimate
$\widehat{\bm{\xi}}\in\mathbb{R}^{NM\times 1}$ obtained at the output of
Algorithm 1. Intuitively, this can lead to unambiguous estimates of $\tau$
since no ambiguity exists in PN covariance $\mathbf{R}(\tau)$ with respect to
$\tau$.

###  G-A PN Observation Model and ML Estimator of Delay

The observation model for the above estimation problem can be written as

$\displaystyle\widehat{\bm{\xi}}=\bm{\xi}+\mathbf{n}\leavevmode\nobreak\ ,$
(119)

where $\bm{\xi}\in\mathbb{R}^{NM\times 1}$ denotes the true PN vector with
$\bm{\xi}\sim{\mathcal{N}}({{\bm{0}}},\mathbf{R}(\tau))$ and
$\mathbf{n}\in\mathbb{R}^{NM\times 1}$ is the estimation noise, independent of
$\bm{\xi}$, whose statistics are given by
$\mathbf{n}\sim{\mathcal{N}}({{\bm{0}}},\mathbf{\Sigma}(\tau))$. Here,
$\mathbf{\Sigma}(\tau)\in\mathbb{R}^{NM\times NM}$ represents the covariance
matrix of the PN estimation error in Algorithm 1 and can be set to the CRB
matrix on PN estimation for the hybrid ML/MAP problem (38). Accordingly, the
statistics of $\widehat{\bm{\xi}}$ in (119) are

$\displaystyle\widehat{\bm{\xi}}\sim{\mathcal{N}}({{\bm{0}}},\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau))\leavevmode\nobreak\
.$ (120)

The ML estimator of $\tau$ from the PN estimates in (120) can be obtained as
[70]

$\displaystyle{\widehat{\tau}}=\arg\min_{\tau}$
$\displaystyle\leavevmode\nobreak\
\Big{\\{}\widehat{\bm{\xi}}^{T}\big{(}\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)\big{)}^{-1}\widehat{\bm{\xi}}$
(121) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\log\det\big{(}\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)\big{)}\Big{\\}}\leavevmode\nobreak\
.$

We identify several challenges in solving (121). First, the inverse of the
covariance matrix of $\widehat{\bm{\xi}}$ needs to be calculated, which leads
to a high computational burden. Second, only a single sample of PN estimate is
available to solve (121), i.e., the sample covariance matrix
$\widehat{\bm{\xi}}\widehat{\bm{\xi}}^{T}$ is rank-one, leading to poor
estimates of $\tau$. Additionally, the formulation in (121) does not allow us
to exploit redundancy in $\mathbf{R}(\tau)$ stemming from its Toeplitz-block
Toeplitz structure in (23).

###  G-B Covariance Matching Approach

To overcome the aforementioned challenges, we propose a covariance matching
(i.e., parametric covariance matrix reconstruction) approach that fits the
sample covariance to the analytical model using a weighted least squares (WLS)
formulation [66]:

$\displaystyle{\widehat{\tau}}=\arg\min_{\tau}$
$\displaystyle\leavevmode\nobreak\
\Big{\\{}\big{[}{\rm{vec}}\big{(}\widehat{\mathbf{R}}\big{)}-{\rm{vec}}\big{(}\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)\big{)}\big{]}^{T}\big{(}\widehat{\mathbf{R}}^{T}\otimes\widehat{\mathbf{R}}\big{)}^{-1}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\times\big{[}{\rm{vec}}\big{(}\widehat{\mathbf{R}}\big{)}-{\rm{vec}}\big{(}\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)\big{)}\big{]}\Big{\\}}\leavevmode\nobreak\
,$ (122)

where $\widehat{\mathbf{R}}=\widehat{\bm{\xi}}\widehat{\bm{\xi}}^{T}$ is the
sample covariance matrix. By the extended invariance principle (EXIP), the
proposed covariance matching based estimator (122) is asymptotically
equivalent to the original ML estimator (121) at high SNRs [66]. Since
$\widehat{\mathbf{R}}$ is rank-one,
$\big{(}\widehat{\mathbf{R}}^{T}\otimes\widehat{\mathbf{R}}\big{)}$ is also
rank-one and thus not invertible. Hence, we resort to the LS reformulation of
(122):

$\displaystyle{\widehat{\tau}}=\arg\min_{\tau}$
$\displaystyle\leavevmode\nobreak\
\left\lVert{\rm{vec}}\big{(}\widehat{\mathbf{R}}\big{)}-{\rm{vec}}\big{(}\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)\big{)}\right\rVert^{2}\leavevmode\nobreak\
.$ (123)

###  G-C High-SNR Approximation

The Toeplitz-block Toeplitz structure of $\mathbf{R}(\tau)$ cannot be
exploited in (123) since $\mathbf{\Sigma}(\tau)$ can have a generic covariance
structure. Moreover, the sample covariance $\widehat{\mathbf{R}}$ is
constructed using a single observation $\widehat{\bm{\xi}}$ and thus
constitutes a very inaccurate estimate of the true covariance
$\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)$. To deal with these issues, we
propose to make the approximation

$\displaystyle\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)\approx\mathbf{R}(\tau)\leavevmode\nobreak\
,$ (124)

which is valid at high SNRs since the CRB matrix elements are inversely
proportional to SNR (i.e., PN estimation becomes sufficiently accurate at high
SNRs such that $\mathbf{R}(\tau)+\mathbf{\Sigma}(\tau)$ is dominated by the PN
covariance $\mathbf{R}(\tau)$). By virtue of this high-SNR approximation, the
problem (123) becomes

$\displaystyle{\widehat{\tau}}=\arg\min_{\tau}$
$\displaystyle\leavevmode\nobreak\
\left\lVert{\rm{vec}}\big{(}\widehat{\mathbf{R}}\big{)}-{\rm{vec}}\big{(}\mathbf{R}(\tau)\big{)}\right\rVert^{2}\leavevmode\nobreak\
,$ (125)

which coincides with the formulation in (56). As seen from (57) and (58), the
Toeplitz-block Toeplitz structure of $\mathbf{R}(\tau)$ can now be effectively
exploited in (125) to estimate delay from a single observation of the PN
vector.

## References

  * [1] K. V. Mishra _et al._ , “Toward millimeter-wave joint radar communications: A signal processing perspective,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 100–114, Sep. 2019.
  * [2] F. Liu _et al._ , “Joint radar and communication design: Applications, state-of-the-art, and the road ahead,” _IEEE Transactions on Communications_ , vol. 68, no. 6, pp. 3834–3862, 2020.
  * [3] T. Wild _et al._ , “Joint design of communication and sensing for beyond 5G and 6G systems,” _IEEE Access_ , vol. 9, pp. 30 845–30 857, 2021\.
  * [4] D. Ma _et al._ , “Joint radar-communication strategies for autonomous vehicles: Combining two key automotive technologies,” _IEEE Signal Processing Magazine_ , vol. 37, no. 4, pp. 85–97, 2020.
  * [5] J. A. Zhang _et al._ , “Enabling joint communication and radar sensing in mobile networks—a survey,” _IEEE Communications Surveys Tutorials_ , vol. 24, no. 1, pp. 306–345, 2022.
  * [6] F. Liu _et al._ , “Integrated sensing and communications: Toward dual-functional wireless networks for 6G and beyond,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 6, pp. 1728–1767, 2022.
  * [7] L. Zheng _et al._ , “Radar and communication coexistence: An overview: A review of recent methods,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 85–99, 2019.
  * [8] A. Hassanien _et al._ , “Dual-function radar communication systems: A solution to the spectrum congestion problem,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 115–126, Sep. 2019.
  * [9] F. Liu _et al._ , “Toward dual-functional radar-communication systems: Optimal waveform design,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 16, pp. 4264–4279, Aug 2018.
  * [10] J. Qian _et al._ , “Radar and communication spectral coexistence in range-dependent interference,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 5891–5906, 2021.
  * [11] W. Saad _et al._ , “A vision of 6G wireless systems: Applications, trends, technologies, and open research problems,” _IEEE Network_ , vol. 34, no. 3, pp. 134–142, 2020.
  * [12] C. Sturm _et al._ , “Waveform design and signal processing aspects for fusion of wireless communications and radar sensing,” _Proceedings of the IEEE_ , vol. 99, no. 7, pp. 1236–1259, July 2011.
  * [13] M. Bică _et al._ , “Generalized multicarrier radar: Models and performance,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 17, pp. 4389–4402, Sep. 2016.
  * [14] F. Zhang _et al._ , “Joint range and velocity estimation with intrapulse and intersubcarrier Doppler effects for OFDM-based RadCom systems,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 662–675, 2020.
  * [15] M. F. Keskin _et al._ , “Limited feedforward waveform design for OFDM dual-functional radar-communications,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 2955–2970, 2021.
  * [16] A. Leshem _et al._ , “Phase noise compensation for OFDM systems,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 21, pp. 5675–5686, 2017.
  * [17] H. Wymeersch _et al._ , “Integration of communication and sensing in 6G: a joint industrial and academic perspective,” in _2021 IEEE 32nd Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC)_ , 2021, pp. 1–7.
  * [18] F. Bozorgi _et al._ , “RF front-end challenges for joint communication and radar sensing,” in _1st IEEE Int. Online Symp. Joint Commun. Sens._ , Feb. 2021.
  * [19] O. H. Salim _et al._ , “Channel, phase noise, and frequency offset in OFDM systems: Joint estimation, data detection, and hybrid cramér-rao lower bound,” _IEEE Transactions on Communications_ , vol. 62, no. 9, pp. 3311–3325, 2014.
  * [20] M. F. Keskin _et al._ , “MIMO-OFDM joint radar-communications: Is ICI friend or foe?” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, no. 6, pp. 1393–1408, 2021.
  * [21] M. Chung _et al._ , “Phase-noise compensation for OFDM systems exploiting coherence bandwidth: Modeling, algorithms, and analysis,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2021.
  * [22] D. D. Lin _et al._ , “The variational inference approach to joint data detection and phase noise estimation in OFDM,” _IEEE Transactions on Signal Processing_ , vol. 55, no. 5, pp. 1862–1874, 2007.
  * [23] R. Wang _et al._ , “Channel estimation, carrier recovery, and data detection in the presence of phase noise in OFDM relay systems,” _IEEE Transactions on Wireless Communications_ , vol. 15, no. 2, pp. 1186–1205, 2016.
  * [24] J. Rodríguez-Fernández, “Joint synchronization and compressive channel estimation for frequency-selective hybrid mmwave MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 1, pp. 548–562, 2022.
  * [25] D. Petrovic _et al._ , “Effects of phase noise on OFDM systems with and without PLL: Characterization and compensation,” _IEEE Transactions on Communications_ , vol. 55, no. 8, pp. 1607–1616, 2007.
  * [26] T. Schenk _et al._ , “On the influence of phase noise induced ICI in MIMO OFDM systems,” _IEEE Communications Letters_ , vol. 9, no. 8, pp. 682–684, 2005.
  * [27] P. Mathecken _et al._ , “Performance analysis of OFDM with Wiener phase noise and frequency selective fading channel,” _IEEE Transactions on Communications_ , vol. 59, no. 5, pp. 1321–1331, 2011.
  * [28] S. Wu _et al._ , “A phase noise suppression algorithm for OFDM-based WLANs,” _IEEE Communications Letters_ , vol. 6, no. 12, pp. 535–537, 2002\.
  * [29] D. D. Lin _et al._ , “Joint estimation of channel response, frequency offset, and phase noise in OFDM,” _IEEE Transactions on Signal Processing_ , vol. 54, no. 9, pp. 3542–3554, 2006.
  * [30] Q. Zou _et al._ , “Compensation of phase noise in OFDM wireless systems,” _IEEE Transactions on Signal Processing_ , vol. 55, no. 11, pp. 5407–5424, 2007.
  * [31] Z. Xue _et al._ , “OFDM radar and communication joint system using opto-electronic oscillator with phase noise degradation analysis and mitigation,” _Journal of Lightwave Technology_ , pp. 1–1, 2022.
  * [32] B. Schweizer _et al._ , “On hardware implementations of stepped-carrier OFDM radars,” in _2018 IEEE/MTT-S International Microwave Symposium - IMS_ , 2018, pp. 891–894.
  * [33] Z. Wei _et al._ , “Toward multi-functional 6G wireless networks: Integrating sensing, communication, and security,” _IEEE Communications Magazine_ , vol. 60, no. 4, pp. 65–71, 2022.
  * [34] A. Tang _et al._ , “Self-interference-resistant IEEE 802.11ad-based joint communication and automotive radar design,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, no. 6, pp. 1484–1499, 2021.
  * [35] Y. Zhuo _et al._ , “Multi-beam joint communication and radar sensing: Beamforming design and interference cancellation,” _IEEE Communications Letters_ , 2022.
  * [36] P. Kumari _et al._ , “Adaptive and fast combined waveform-beamforming design for mmwave automotive joint communication-radar,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, no. 4, pp. 996–1012, 2021\.
  * [37] M. A. Uusitalo _et al._ , “6G vision, value, use cases and technologies from european 6G flagship project Hexa-X,” _IEEE Access_ , vol. 9, pp. 160 004–160 020, 2021.
  * [38] R. F. Tigrek _et al._ , “OFDM signals as the radar waveform to solve Doppler ambiguity,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 48, no. 1, pp. 130–143, Jan 2012.
  * [39] M. Braun, “OFDM radar algorithms in mobile communication networks,” _Karlsruher Institutes für Technologie_ , 2014.
  * [40] S. Mercier _et al._ , “Comparison of correlation-based OFDM radar receivers,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 56, no. 6, pp. 4796–4813, 2020.
  * [41] R. Xie _et al._ , “Performance analysis of joint range-velocity estimator with 2D-MUSIC in OFDM radar,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 4787–4800, 2021.
  * [42] M. C. Budge _et al._ , “Range correlation effects in radars,” in _The Record of the 1993 IEEE National Radar Conference_ , April 1993, pp. 212–216.
  * [43] A. R. Chiriyath _et al._ , “Joint radar-communications information bounds with clutter: The phase noise menace,” in _2016 IEEE Radar Conference (RadarConf)_ , May 2016, pp. 1–6.
  * [44] M. Gerstmair _et al._ , “On the safe road toward autonomous driving: Phase noise monitoring in radar sensors for functional safety compliance,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 60–70, 2019.
  * [45] C. Aydogdu _et al._ , “Radar interference mitigation for automated driving: Exploring proactive strategies,” _IEEE Signal Processing Magazine_ , vol. 37, no. 4, pp. 72–84, 2020.
  * [46] S. Wu _et al._ , “OFDM systems in the presence of phase noise: consequences and solutions,” _IEEE Transactions on Communications_ , vol. 52, no. 11, pp. 1988–1996, 2004.
  * [47] A. A. Nasir _et al._ , “Phase noise in MIMO systems: Bayesian Cramér–Rao bounds and soft-input estimation,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 10, pp. 2675–2692, 2013.
  * [48] A. Demir, “Computing timing jitter from phase noise spectra for oscillators and phase-locked loops with white and $1/f$ noise,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 53, no. 9, pp. 1869–1884, Sep. 2006.
  * [49] A. Chorti _et al._ , “A spectral model for RF oscillators with power-law phase noise,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 53, no. 9, pp. 1989–1999, Sep. 2006.
  * [50] J. Tao _et al._ , “Estimation of channel transfer function and carrier frequency offset for OFDM systems with phase noise,” _IEEE Transactions on Vehicular Technology_ , vol. 58, no. 8, pp. 4380–4387, 2009.
  * [51] G. Hakobyan _et al._ , “A novel intercarrier-interference free signal processing scheme for OFDM radar,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 6, pp. 5158–5167, 2017.
  * [52] S. Bay _et al._ , “On the hybrid Cramér Rao bound and its application to dynamical phase estimation,” _IEEE Signal Processing Letters_ , vol. 15, pp. 453–456, 2008.
  * [53] Y. Noam _et al._ , “Notes on the tightness of the hybrid Cramér–Rao lower bound,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 6, pp. 2074–2084, 2009.
  * [54] Z. Zhang _et al._ , “Full duplex techniques for 5G networks: self-interference cancellation, protocol design, and relay selection,” _IEEE Communications Magazine_ , vol. 53, no. 5, pp. 128–137, 2015.
  * [55] P. Kumari _et al._ , “IEEE 802.11ad-based radar: An approach to joint vehicular communication-radar system,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 4, pp. 3012–3027, April 2018.
  * [56] C. B. Barneto _et al._ , “Full-duplex OFDM radar with LTE and 5G NR waveforms: Challenges, solutions, and measurements,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 67, no. 10, pp. 4042–4054, 2019.
  * [57] A. Demir _et al._ , “Phase noise in oscillators: a unifying theory and numerical methods for characterization,” _IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications_ , vol. 47, no. 5, pp. 655–674, 2000.
  * [58] P. Rabiei _et al._ , “A non-iterative technique for phase noise ICI mitigation in packet-based OFDM systems,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 11, pp. 5945–5950, 2010.
  * [59] X. Quan _et al._ , “Impacts of phase noise on digital self-interference cancellation in full-duplex communications,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 7, pp. 1881–1893, 2017.
  * [60] C. R. Berger _et al._ , “Signal processing for passive radar using OFDM waveforms,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 4, no. 1, pp. 226–238, 2010.
  * [61] L. Zheng _et al._ , “Super-resolution delay-Doppler estimation for OFDM passive radar,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 9, pp. 2197–2210, May 2017.
  * [62] L. Gaudio _et al._ , “On the effectiveness of OTFS for joint radar parameter estimation and communication,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 9, pp. 5951–5965, 2020.
  * [63] A. Murat _et al._ , “Phase-noise-induced performance limits for DPSK modulation with and without frequency feedback,” _Journal of Lightwave Technology_ , vol. 11, no. 2, pp. 290–302, 1993.
  * [64] M. Wax _et al._ , “Efficient inversion of Toeplitz-block Toeplitz matrix,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 31, no. 5, pp. 1218–1221, 1983.
  * [65] P. Stoica _et al._ , “On reparametrization of loss functions used in estimation and the invariance principle,” _Signal Processing_ , vol. 17, no. 4, pp. 383–387, 1989.
  * [66] B. Ottersten _et al._ , “Covariance matching estimation techniques for array signal processing applications,” _Digital Signal Processing_ , vol. 8, no. 3, pp. 185–210, 1998.
  * [67] X. Du _et al._ , “Toeplitz structured covariance matrix estimation for radar applications,” _IEEE Signal Processing Letters_ , vol. 27, pp. 595–599, 2020.
  * [68] A. Aubry _et al._ , “A new sequential optimization procedure and its applications to resource allocation for wireless systems,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 24, pp. 6518–6533, 2018.
  * [69] T. F. Chan, “An optimal circulant preconditioner for Toeplitz systems,” _SIAM Journal on Scientific and Statistical Computing_ , vol. 9, no. 4, pp. 766–771, 1988.
  * [70] A. Swindlehurst _et al._ , “Maximum likelihood methods in radar array signal processing,” _Proceedings of the IEEE_ , vol. 86, no. 2, pp. 421–441, 1998.