Position and Orientation Estimation through Millimeter Wave MIMO in 5G Systems

TL;DR

This paper investigates the potential of millimeter wave MIMO systems in 5G for accurate position and orientation estimation, deriving theoretical bounds and proposing a novel estimation algorithm that performs well across various conditions.

Contribution

It introduces a new two-stage estimation algorithm for position and rotation angle that approaches the theoretical CRB in millimeter wave 5G systems with a single transmitter.

Findings

Accurate position and orientation estimation is feasible with a single transmitter.

The proposed algorithm achieves near-CRB performance at high SNR.

Effective in line-of-sight and obstructed conditions.

Abstract

Millimeter wave signals and large antenna arrays are considered enabling technologies for future 5G networks. While their benefits for achieving high-data rate communications are well-known, their potential advantages for accurate positioning are largely undiscovered. We derive the Cram\'{e}r-Rao bound (CRB) on position and rotation angle estimation uncertainty from millimeter wave signals from a single transmitter, in the presence of scatterers. We also present a novel two-stage algorithm for position and rotation angle estimation that attains the CRB for average to high signal-to-noise ratio. The algorithm is based on multiple measurement vectors matching pursuit for coarse estimation, followed by a refinement stage based on the space-alternating generalized expectation maximization algorithm. We find that accurate position and rotation angle estimation is possible using signals from a single transmitter, in either line-of- sight, non-line-of-sight, or obstructed-line-of-sight conditions.

Tables

Table 1. TABLE I: Unknown conditions

SNR (in dB)	-20	-10	0	10
$\Delta v$	5.5	5.2	5	5.3

Equations

$$\mathbf{H}[n]={\mathbf{A}}_{\mathrm{Rx}}[n]\mathbf{\Gamma }[n]{\mathbf{A}}_{\mathrm{Tx}}^{\mathrm{H}}[n],$$

$${\mathbf{A}}_{\mathrm{Tx}}[n]$$

$${\mathbf{A}}_{\mathrm{Rx}}[n]$$

$$\mathbf{\Gamma }[n]=\sqrt{N_t{N}_r}\mathrm{diag}\left\{\frac{h_0}{\sqrt{{\rho }_0}}{e}^{-j2\pi n{\tau }_0/(N{T}_s)},\dots ,\frac{h_K}{\sqrt{{\rho }_K}}{e}^{-j2\pi n{\tau }_K/(N{T}_s)}\right\},$$

$${\mathbf{a}}_{\mathrm{Tx},n}({\theta }_{\mathrm{Tx},k})=$$

$$\frac{1}{\sqrt{N_t}}[e^{-j\frac{N_t-1}{2}\frac{2\pi }{{\lambda }_n}d\sin ({\theta }_{\mathrm{Tx},k})},\dots ,e^{j\frac{N_t-1}{2}\frac{2\pi }{{\lambda }_n}d\sin ({\theta }_{\mathrm{Tx},k})}{]}^{\mathrm{T}},$$

$${\mathbf{y}}^{(g)}[n]=\mathbf{H}[n]{\mathbf{F}}^{(g)}[n]{\mathbf{x}}^{(g)}[n]+{\mathbf{n}}^{(g)}[n],$$

$$\mathbf{\eta }={\left[\begin{array}{c}{\mathbf{\eta }}_0^{\mathrm{T}},\dots ,{\mathbf{\eta }}_K^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}},$$

$${\mathbf{\eta }}_k={\left[\begin{array}{c}{\tau }_k,{\mathbf{\theta }}_k^{\mathrm{T}},{\tilde{\mathbf{h}}}_k^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}},$$

$${\mathbb{E}}_{\mathbf{y}|\mathbf{\eta }}\left[(\hat{\mathbf{\eta }}-\mathbf{\eta })(\hat{\mathbf{\eta }}-\mathbf{\eta }{)}^{\mathrm{T}}\right]⪰{\mathbf{J}}_{\mathbf{\eta }}^{-1},$$

$${\mathbf{J}}_{\mathbf{\eta }}\triangleq {\mathbb{E}}_{\mathbf{y}|\mathbf{\eta }}\left[-\frac{{\partial }^2\ln f(\mathbf{y}|\mathbf{\eta })}{\partial \mathbf{\eta }\partial {\mathbf{\eta }}^T}\right],$$

$$f(\mathbf{y}|\mathbf{\eta })\propto \exp \left\{\frac{2}{N_0}\sum _{n=0}^{N-1}\Re \{{\mathbf{\mu }}^{\mathrm{H}}[n]\mathbf{y}[n]\}-\frac{1}{N_0}\sum _{n=0}^{N-1}\Vert \mathbf{\mu }[n]{\Vert }_2^2\right\},$$

$${\mathbf{J}}_{\mathbf{\eta }}=\left[\begin{array}{ccc}\mathbf{\Psi }({\mathbf{\eta }}_0,{\mathbf{\eta }}_0)& \dots & \mathbf{\Psi }({\mathbf{\eta }}_0,{\mathbf{\eta }}_K)\\ ⋮& \ddots & ⋮\\ \mathbf{\Psi }({\mathbf{\eta }}_K,{\mathbf{\eta }}_0)& \dots & \mathbf{\Psi }({\mathbf{\eta }}_K,{\mathbf{\eta }}_K)\end{array}\right],$$

$$\mathbf{\Psi }({\mathbf{x}}_r,{\mathbf{x}}_s)\triangleq {\mathbb{E}}_{\mathbf{y}|\mathbf{\eta }}\left[-\frac{{\partial }^2\ln f(\mathbf{y}|\mathbf{\eta })}{\partial {\mathbf{x}}_r\partial {\mathbf{x}}_s^{\mathrm{T}}}\right].$$

$$\mathbf{\Psi }({\mathbf{\eta }}_r,{\mathbf{\eta }}_s)=\left[\begin{array}{ccc}\Psi ({\tau }_r,{\tau }_s)& \mathbf{\Psi }({\tau }_r,{\mathbf{\theta }}_s)& \mathbf{\Psi }({\tau }_r,{\mathbf{h}}_s)\\ \mathbf{\Psi }({\mathbf{\theta }}_r,{\tau }_s)& \mathbf{\Psi }({\mathbf{\theta }}_r,{\mathbf{\theta }}_s)& \mathbf{\Psi }({\mathbf{\theta }}_r,{\mathbf{h}}_s)\\ \mathbf{\Psi }({\mathbf{h}}_r,{\tau }_s)& \mathbf{\Psi }({\mathbf{h}}_r,{\mathbf{\theta }}_s)& \mathbf{\Psi }({\mathbf{h}}_r,{\mathbf{h}}_s)\end{array}\right].$$

$${\mathbf{J}}_{\tilde{\mathbf{\eta }}}=\mathbf{T}{\mathbf{J}}_{\mathbf{\eta }}{\mathbf{T}}^{\mathrm{T}},$$

$$\mathbf{T}\triangleq \frac{\partial {\mathbf{\eta }}^{\mathrm{T}}}{\partial \tilde{\mathbf{\eta }}}.$$

$${\tau }_0$$

$${\tau }_k$$

$${\theta }_{\mathrm{Tx},0}$$

$${\theta }_{\mathrm{Tx},k}$$

$${\theta }_{\mathrm{Rx},k}$$

$${\theta }_{\mathrm{Rx},0}$$

$$\mathbf{T}=\left[\begin{array}{ccc}{\mathbf{T}}_{0,0}& \dots & {\mathbf{T}}_{K,0}\\ ⋮& \ddots & ⋮\\ {\mathbf{T}}_{0,K}& \dots & {\mathbf{T}}_{K,K}\end{array}\right],$$

$${\mathbf{T}}_{k,k^{\prime }}\triangleq \frac{\partial {\mathbf{\eta }}_k^{\mathrm{T}}}{\partial {\tilde{\mathbf{\eta }}}_{k^{\prime }}}.$$

$${\mathbf{T}}_{k,k^{\prime }}=\left[\begin{array}{ccc}\partial {\tau }_k/\partial {\mathbf{s}}_{k^{\prime }}& \partial {\mathbf{\theta }}_k^{\mathrm{T}}/\partial {\mathbf{s}}_{k^{\prime }}& \partial {\tilde{\mathbf{h}}}_k^{\mathrm{T}}/\partial {\mathbf{s}}_{k^{\prime }}\\ \partial {\tau }_k/\partial {\tilde{\mathbf{h}}}_{k^{\prime }}& \partial {\mathbf{\theta }}_k^{\mathrm{T}}/\partial {\tilde{\mathbf{h}}}_{k^{\prime }}& \partial {\tilde{\mathbf{h}}}_k^{\mathrm{T}}/\partial {\tilde{\mathbf{h}}}_{k^{\prime }}\end{array}\right],$$

$${\mathbf{T}}_{k,0}=\left[\begin{array}{ccc}\partial {\tau }_k/\partial \mathbf{p}& \partial {\mathbf{\theta }}_k^{\mathrm{T}}/\partial \mathbf{p}& \partial {\tilde{\mathbf{h}}}_k^{\mathrm{T}}/\partial \mathbf{p}\\ \partial {\tau }_k/\partial \alpha & \partial {\mathbf{\theta }}_k^{\mathrm{T}}/\partial \alpha & \partial {\tilde{\mathbf{h}}}_k^{\mathrm{T}}/\partial \alpha \\ \partial {\tau }_k/\partial {\tilde{\mathbf{h}}}_0& \partial {\mathbf{\theta }}_k^{\mathrm{T}}/\partial {\tilde{\mathbf{h}}}_0& \partial {\tilde{\mathbf{h}}}_k^{\mathrm{T}}/\partial {\tilde{\mathbf{h}}}_0\end{array}\right],$$

$$\partial {\tau }_0/\partial \mathbf{p}$$

$$\partial {\theta }_{\mathrm{Tx},0}/\partial \mathbf{p}$$

$$\partial {\theta }_{\mathrm{Rx},0}/\partial \mathbf{p}$$

$$\partial {\theta }_{\mathrm{Rx},k}/\partial \alpha$$

$$\partial {\tau }_k/\partial \mathbf{p}$$

$$\partial {\tau }_k/\partial {\mathbf{s}}_k$$

$$\partial {\theta }_{\mathrm{Tx},k}/\partial {\mathbf{s}}_k$$

$$\partial {\theta }_{\mathrm{Rx},k}/\partial \mathbf{p}$$

$$\partial {\theta }_{\mathrm{Rx},k}/\partial {\mathbf{s}}_k$$

$$\mathrm{PEB}=\sqrt{\mathrm{tr}\left\{[{\mathbf{J}}_{\tilde{\mathbf{\eta }}}^{-1}{]}_{1:2,1:2}\right\}},$$

Full text

I Introduction

Fifth generation (5G) communication networks will likely adopt millimeter-wave (mm-wave) and massive multiple-input-multiple-output (MIMO) technologies, thanks to a number of favorable properties. In particular, operating at carrier frequencies beyond 30 GHz, with large available bandwidths, mm-wave can provide extremely high data rates to users through dense spatial multiplexing by using a large number of antennas [2, 3]. While these properties are desirable for 5G services, mm-wave communications also face a number of challenges. Among these, the severe path loss at high carrier frequencies stands out. The resulting loss in signal-to-noise ratio (SNR) must be compensated through sophisticated beamforming at the transmitter and/or receiver side, leading to highly directional links [4, 5, 6]. However, beamforming requires knowledge of the propagation channel. Significant progress has been made in mm-wave channel estimation, by exploiting sparsity and related compressed sensing tools, such as distributed compressed sensing-simultaneous orthogonal matching pursuit (DCS-SOMP) [7], compressive sampling matched pursuit (CoSOMP) [8], and group sparse compressed sensing (GCS) [9]. In particular, since at mm-wave frequencies only the line-of-sight (LOS) path and a few dominant multi-path components contribute to the received power, mm-wave channels are sparse in the angular domain [10, 11]. This is because in mm-wave frequencies, the received power of diffuse scattering and multiple-bounce specular reflections are much lower than that in LOS and single-bounce specular reflection [12, 13, 14]. Different compressed sensing (CS) methods for mm-wave channel estimation are proposed in[15, 16, 17, 18, 19, 20, 21, 22, 23]. In [15], a method for the estimation of angle-of-arrival (AOA), angle-of-departure (AOD), and channel gains is proposed based on the compressive beacons on the downlink. A method for the continuous estimation of mm-wave channel parameters is proposed in [22], while [21] applies CS tools with refinement in the angular domain. In [16], a CS method is proposed based on the redundant dictionary matrices. A two-stage algorithm with one-time feedback that is robust to noise is used in [20]. In [17], an adaptive CS method is proposed based on a hierarchical multi-resolution codebook design for the estimation of single-path and multi-path mm-wave channels. In [18], a beam selection procedure for the multiuser mm-wave MIMO channels with analog beamformers is proposed. In [19], a CS approach with reduced training overhead was considered. Finally, CS tools are used in [23] for the sparse estimation of power angle profiles of the mm-wave channels and compared with the codebook designs in terms of overhead reduction. However, in all the aforementioned papers, a narrow-band mm-wave channel model is used. When the bandwidths becomes larger, one needs to consider the effect of the delays of different paths in the mm-wave channel model, i.e., the wide-band mm-wave channel model.

Channel estimation provides information of the AOA/AOD and thus of the relative location of the transmitter and receiver. In addition, location information can serve as a proxy for channel information to perform beamforming: when the location of the user is known, the base station (BS) can steer its transmission to the user, either directly or through a reflected path. This leads to synergies between localization and communication. The use of 5G technologies to obtain position and orientation was previously explored in [24, 25, 26] for mm-wave and in [27, 28, 29] for massive MIMO. The early work [24] considered estimation and tracking of AOA through beam-switching. User localization was treated in [25], formulated as a hypothesis testing problem, limiting the spatial resolution. A different approach was taken in [26], where meter-level positioning accuracy was obtained by measuring received signal strength levels. A location-aided beamforming method was proposed in [30] to speed up initial access between nodes. In the massive MIMO case, [27] considered the estimation of angles, [31] proposed a direct localization method by jointly processing the observations at the distributed massive MIMO BSs, while [28] treated the joint estimation of delay, AOD, and AOA, in the LOS conditions and evaluated the impact of errors in delays and phase shifters, and [1] derived sufficient conditions for a nonsingular Fisher information matrix (FIM) of delay, AOD, AOA, and channel coefficients. A hybrid time-difference-of-arrival (TDOA), AOA, and AOD localization was proposed in [32] using linearization. In [29], positioning was solved using a Gaussian process regressor, operating on a vector of received signal strengths through fingerprinting. While latter this approach is able to exploit non-line-of-sight (NLOS) propagation, it does not directly harness the geometry of the environment. Complementarily to the use of mm-wave frequencies, approaches for localization using centimeter-wave (cm-wave) signals have been recently proposed as well. The combination of TDOAs and AOAs using an extended Kalman filter (EKF) was presented in [33, 34], where the mobile station (MS) has a single antenna, while the BS employs an antenna array. This method assumes LOS propagation thanks to the high density of access nodes and provides sub-meter accuracy even for moving devices.

In this paper, we show that mm-wave and large MIMO are enabling technologies for accurate positioning and device orientation estimation with only one BS, even when the LOS path is blocked. The limited scattering and high-directivity are unique characteristics of the mm-wave channel and large MIMO systems, respectively. We derive fundamental bounds on the position and orientation estimation accuracy, for LOS111LOS is defined as the condition where the LOS path exists and there are no scatterers., NLOS222NLOS is defined as the condition where there are scatterers and the LOS path is not blocked., and obstructed-line-of-sight (OLOS)333OLOS is referred to the condition where the LOS path is blocked and only the signals from the scatterers are received. conditions. These bounds indicate that the information from the NLOS links help to estimate the location and orientation of the MS. We also propose a novel three-stage position and orientation estimation technique, which is able to attain the bounds at average to high SNR. The first stage of the technique harnesses sparsity of the mm-wave channel in the AOA and AOD domain [10, 11]. Moreover, the sparsity support does not vary significantly with frequency, allowing us to use DCS-SOMP across different carriers. The delay can then be estimated on a per-path basis. As DCS-SOMP limits the AOA and AOD to a predefined grid, we propose a refinement stage, based on the space-alternating generalized expectation maximization (SAGE) algorithm. Finally, in the last stage, we employ a least-squares approach with extended invariance principle (EXIP) to recover position and orientation [35, 36].

II System Model

We consider a MIMO system with a BS equipped with $N_{t}$ antennas and a MS equipped by $N_{r}$ antennas operating at a carrier frequency $f_{c}$ (corresponding to wavelength $\lambda_{c}$) and bandwidth $B$. Locations of the MS and BS are denoted by $\mathbf{p}=[p_{x},p_{y}]^{\mathrm{T}}\in\mathbb{R}^{2}$ and $\mathbf{q}=[q_{x},q_{y}]^{\mathrm{T}}\in\mathbb{R}^{2}$ with the $\alpha\in[0,2\pi)$ denoting the rotation angle of the MS’s antenna array. The value of $\mathbf{q}$ is assumed to be known, while $\mathbf{p}$ and $\alpha$ are unknown.

II-A Transmitter Model

We consider the transmission of orthogonal frequency division multiplexing (OFDM) signals as in [37], where a BS with hybrid analog/digital precoder communicates with a single MS. At the BS, $G$ signals are transmitted sequentially, where the $g$-th transmission comprises $M_{t}$ simultaneously transmitted symbols $\mathbf{x}^{(g)}[n]=[x_{1}[n],\ldots,x_{M_{t}}[n]]^{\mathrm{T}}\in\mathbb{C}^{M_{t}}$ for each subcarrier $n=0,\ldots,N-1$. The symbols are first precoded and then transformed to the time-domain using $N$-point inverse fast Fourier transform (IFFT). A cyclic prefix (CP) of length $T_{\mathrm{CP}}=DT_{s}$ is added before applying the radio-frequency (RF) precoding where $D$ is the length of CP in symbols. Here, $T_{s}=1/B$ denotes the sampling period and $T_{\mathrm{CP}}$ is assumed to exceed the delay spread of the channel. The transmitted signal over subcarrier $n$ at time $g$ can be expressed as $\mathbf{F}^{(g)}[n]\mathbf{x}^{(g)}[n]$. The beamforming matrix $\mathbf{F}[n]\in\mathbb{C}^{N_{t}\times M_{t}}$ is defined as $\mathbf{F}[n]=\mathbf{F}_{\mathrm{RF}}\mathbf{F}_{\mathrm{BB}}[n]$ where $\mathbf{F}_{\mathrm{RF}}$ is implemented using the analog phase shifters with the entries of the form $e^{j\phi_{m,n}}$, where $\{\phi_{m,n}\}$ are given phases, and $\mathbf{F}_{\mathrm{BB}}[n]$ is the digital beamformer, and overall they satisfy a total power constraint $\|\mathbf{F}_{\mathrm{RF}}\mathbf{F}_{\mathrm{BB}}[n]\|_{\mathrm{F}}=1$. Considering the sparsity of the mm-wave channels one usually needs much less beams $M_{t}$ than antenna elements $N_{t}$, i.e., $M_{t}\ll N_{t}$. Also, the presence of $\mathbf{F}[n]$ in the proposed model leads to the extension of system model to multi-user mm-wave downlink systems with a limited feedback channel from MSs to the BS. Our work does not assume any specific beamformer. We will provide general expressions that permit the study of the impact on performance and optimization of different choices of beamformers $\mathbf{F}^{(g)}[n]$ and signals $\mathbf{x}^{(g)}[n]$, although this is out of the scope of the paper. Our approach is also compatible with beam reference signal (initial access) procedures, and it could be complemented with a Bayesian recursive tracker with user-specific precoding.

II-B Channel Model

Fig. 1 shows the position-related parameters of the channel. These parameters include $\theta_{\mathrm{Rx},k}$, $\theta_{\mathrm{Tx},k}$, and $d_{k}=c\tau_{k}$, denoting the AOA, AOD, and the path length (with time-of-arrival (TOA) $\tau_{k}$ and the speed of light $c$) of the $k$-th path ($k=0$ for the LOS path and $k>0$ the NLOS paths). For each NLOS path, there is a scatterer with unknown location $\mathbf{s}_{k}$, for which we define $d_{k,1}=\|\mathbf{s}_{k}-\mathbf{q}\|_{2}$ and $d_{k,2}=\|\mathbf{p}-\mathbf{s}_{k}\|_{2}$. We now introduce the channel model, under a frequency-dependent array response [11], suitable for wideband communication (with fractional bandwidth $B/f_{c}$ up to $50\%$). Assuming $K+1$ paths and a channel that remains constant during the transmission of $G$ symbols, the $N_{r}\times N_{t}$ channel matrix associated with subcarrier $n$ is expressed as

[TABLE]

for response vectors

[TABLE]

and

[TABLE]

for path loss $\rho_{k}$ and complex channel gain $h_{k}$, respectively, of the $k$-th path. For later use, we introduce $\tilde{h}_{k}=\sqrt{(N_{t}N_{r})/\rho_{k}}h_{k}$ and $\gamma_{n}(h_{k},\tau_{k})=\tilde{h}_{k}e^{-j2\pi n\tau_{k}/(NT_{s})}$.

The structure of the frequency-dependent antenna steering and response vectors $\mathbf{a}_{\mathrm{Tx},n}(\theta_{\mathrm{Tx},k})\in\mathbb{C}^{N_{t}}$ and $\mathbf{a}_{\mathrm{Rx},n}(\theta_{\mathrm{Rx},k})\in\mathbb{C}^{N_{r}}$ depends on the specific array structure. For the case of a uniform linear array (ULA), which will be the example studied in this paper, we recall that (the response vector $\mathbf{a}_{\mathrm{Rx},n}(\theta_{\mathrm{Rx},k})$ is obtained similarly)

[TABLE]

where $\lambda_{n}=c/(n/(NT_{s})+f_{c})$ is the signal wavelength at the $n$-th subcarrier and $d$ denotes the distance between the antenna elements (we will use $d=\lambda_{c}/2$). We note that when $B\ll f_{c}$, $\lambda_{n}\approx\lambda_{c}$, and $\eqref{steringvector1}$ reverts to the standard narrow-band model.

II-C Received Signal Model

The received signal for subcarrier $n$ and transmission $g$, after CP removal and fast Fourier transform (FFT), can be expressed as

[TABLE]

where $\mathbf{n}^{(g)}[n]\in\mathbb{C}^{N_{r}}$ is a Gaussian noise vector with zero mean and variance $N_{0}/2$ per real dimension. Our goal is now to estimate the position $\mathbf{p}$ and orientation $\alpha$ of the MS from $\{\mathbf{y}^{(g)}[n]\}_{\forall n,g}$. We will first derive a fundamental lower bound on the estimation uncertainty and then propose a novel practical estimator.

III Position and Orientation Estimation: Fundamental Bounds

In this section, we derive the FIM and the Cramér-Rao bound (CRB) for the estimation problem of position and orientation of the MS for LOS, NLOS, and OLOS. To simplify the notation and without loss of generality, we consider the case of $G=1$, i.e., only 1 OFDM symbol is transmitted.

III-A FIM Derivation for Channel Parameters

Let $\boldsymbol{\eta}\in\mathbb{R}^{5(K+1)}$ be the vector consisting of the unknown channel parameters

[TABLE]

in which $\boldsymbol{\eta}_{k}$ consists of the unknown channel parameters (delay, AOD, AOA, and channel coefficients) for the $k$-th path

[TABLE]

where $\tilde{\mathbf{h}}_{k}=[\tilde{h}_{\mathrm{R},k},\tilde{h}_{\mathrm{I},k}]^{\mathrm{T}}$ contains the real and imaginary parts defined as $\tilde{h}_{\mathrm{R},k}$ and $\tilde{h}_{\mathrm{I},k}$, respectively, and $\boldsymbol{\theta}_{k}=\begin{bmatrix}\theta_{\mathrm{Tx},k},\theta_{\mathrm{Rx},k}\end{bmatrix}^{\mathrm{T}}$.

Defining $\hat{\boldsymbol{\eta}}$ as the unbiased estimator of $\boldsymbol{\eta}$, the mean squared error (MSE) is bounded as [38]

[TABLE]

in which $\mathbb{E}_{\mathbf{y}|\boldsymbol{\eta}}[.]$ denotes the expectation parameterized by the unknown parameters $\boldsymbol{\eta}$, and $\mathbf{J}_{\boldsymbol{\eta}}$ is the $5(K+1)\times 5(K+1)$ FIM defined as

[TABLE]

where $f(\mathbf{y}|\boldsymbol{\eta})$ is the likelihood function of the random vector $\mathbf{y}$ conditioned on $\boldsymbol{\eta}$. More specifically, $f(\mathbf{y}|\boldsymbol{\eta})$ can be written as [39]

[TABLE]

where $\boldsymbol{\mu}[n]\triangleq\mathbf{H}[n]\mathbf{F}[n]\mathbf{x}[n]$ and $\propto$ denotes equality up to irrelevant constants.

The FIM in (10) can be structured as

[TABLE]

in which $\mathbf{\Psi}(\mathbf{x}_{r},\mathbf{x}_{s})$ is defined as

[TABLE]

The $5\times 5$ matrix $\mathbf{\Psi}(\boldsymbol{\eta}_{r},\boldsymbol{\eta}_{s})$ is structured as

[TABLE]

The entries of $\mathbf{\Psi}(\boldsymbol{\eta}_{r},\boldsymbol{\eta}_{s})$ are derived in Appendix A.

III-B FIM for Position and Orientation

We determine the FIM in the position space through a transformation of variables from $\boldsymbol{\eta}$ to $\tilde{\boldsymbol{\eta}}=\begin{bmatrix}\tilde{\boldsymbol{\eta}}^{\mathrm{T}}_{0},\ldots,\tilde{\boldsymbol{\eta}}^{\mathrm{T}}_{K}\end{bmatrix}^{\mathrm{T}}$, where $\tilde{\boldsymbol{\eta}}_{k}=\begin{bmatrix}\mathbf{s}^{\mathrm{T}}_{k},\tilde{\mathbf{h}}^{\mathrm{T}}_{k}\end{bmatrix}^{\mathrm{T}}$ for $k>0$ and $\tilde{\boldsymbol{\eta}}_{0}=\begin{bmatrix}\mathbf{p}^{\mathrm{T}},\alpha,\tilde{\mathbf{h}}^{\mathrm{T}}_{0}\end{bmatrix}^{\mathrm{T}}$. If the LOS path is blocked (i.e., OLOS), we note that we must consider $\boldsymbol{\eta}_{\mathrm{olos}}=[\boldsymbol{\eta}^{\mathrm{T}}_{1},\ldots,\boldsymbol{\eta}^{\mathrm{T}}_{K}]^{\mathrm{T}}$ and $\tilde{\boldsymbol{\eta}}_{\mathrm{olos}}=[\mathbf{p}^{\mathrm{T}},\alpha,\tilde{\boldsymbol{\eta}}^{\mathrm{T}}_{1},\ldots,\tilde{\boldsymbol{\eta}}^{\mathrm{T}}_{K}]^{\mathrm{T}}$.

The FIM of $\tilde{\boldsymbol{\eta}}$ is obtained by means of the $(4K+5)\times 5(K+1)$ transformation matrix $\mathbf{T}$ as

[TABLE]

where

[TABLE]

The entries of $\mathbf{T}$ can be obtained by the relations between the parameters in $\boldsymbol{\eta}$ and $\tilde{\boldsymbol{\eta}}$ from the geometry of the problem shown in Fig. 1 as:

[TABLE]

Consequently, we obtain

[TABLE]

in which $\mathbf{T}_{k,k^{\prime}}$ is defined as

[TABLE]

For $k^{\prime}\neq 0$, $\mathbf{T}_{k,k^{\prime}}$ is obtained as

[TABLE]

and $\mathbf{T}_{k,0}$ is obtained as

[TABLE]

where

[TABLE]

[TABLE]

and $\partial\tilde{\mathbf{h}}^{\mathrm{T}}_{k}/\tilde{\mathbf{h}}_{k}=\mathbf{I}_{2}$ for $k\geq 0$. The rest of entries in $\mathbf{T}$ are zero.

III-C Bounds on Position and Orientation Estimation Error

The position error bound (PEB) is obtained by inverting $\mathbf{J}_{\tilde{\boldsymbol{\eta}}}$, adding the diagonal entries of the $2\times 2$ sub-matrix, and taking the root square as:

[TABLE]

and the rotation error bound (REB) is obtained as:

[TABLE]

where the operations $[.]_{1:2,1:2}$ and $[.]_{3,3}$ denote the selection of the first $2\times 2$ sub-matrix and the third diagonal entry of $\mathbf{J}^{-1}_{\tilde{\boldsymbol{\eta}}}$, respectively.

III-D The Effect of Multi-Path Components on Position and Orientation Estimation Error

In this subsection, we discuss the effect of adding multi-path components (MPCs) for localization under different conditions. As the number of antennas in the MS increases, the scalar product between steering vectors corresponding to different receive directions tends to vanish, i.e. $|\mathbf{a}^{\mathrm{H}}_{\mathrm{Rx},n}(\theta_{\mathrm{Rx},r})\mathbf{a}_{\mathrm{Rx},n}(\theta_{\mathrm{Rx},s})|\ll 1$ for $\theta_{\mathrm{Rx},r}\neq\theta_{\mathrm{Rx},s}$. Also, increasing the number of antenna elements in the transmitter results in narrower beams and the spatial correlation between different beams is reduced. Moreover, as the system bandwidth increases, the different MPCs coming from different scatterers can be more easily resolved. In other words, the MPCs can be considered to be orthogonal [40, 41]. Consequently, large $N_{t}$, $N_{r}$, and bandwidth lead to very small multipath cross-correlation terms in the FIM [42]. Ignoring those terms, the approximate expression for the equivalent Fisher information matrix (EFIM) of position and rotation angle $\mathbf{J}_{e}(\mathbf{p},\alpha)$ with large $N_{t}$, $N_{r}$, and bandwidth is444In computing (29), we used the fact that the last two rows of $\mathbf{T}_{k,0}$ are zero for $k\neq 0$.

[TABLE]

where

[TABLE]

in which $\tilde{\mathbf{T}}_{0,0}$ is the $3\times 3$ sub-matrix in the transformation matrix $\mathbf{T}_{k,0}$ for $k=0$ in (26) containing the derivatives with respect to $\mathbf{p}$ and $\alpha$, $[.]_{1:3,1:3}$ denotes the selection of the first $3\times 3$ sub-matrix, and $\mathbf{\Lambda}_{e,0}$ denotes the EFIM of the delay, AOD, and AOA from LOS, i.e., $\{\tau_{0},\theta_{\mathrm{Tx},0},\theta_{\mathrm{Rx},0}\}$. From simulations, it is observed that the exact and approximate FIM lead to nearly identical PEBs, under the mentioned conditions. Hence, greedy techniques from compressed sensing, which extract path after path, are a natural tool for such scenarios. In the LOS case, (29) only contains the term corresponding to $k=0$, i.e., the first term. When MPCs are present, the terms corresponding to $k\geq 1$ appear, i.e., the second summand in (29), which contains terms that are added and others that are subtracted (because the scatterer location is an additional parameter that has to be estimated for each MPC [43, eq. (3.59)]). The additive terms imply that the presence of MPCs help in the estimation of the MS localization, as they add information to the EFIM. In general the contribution of the MPCs results in a positive contribution to the FIM, and hence in a reduction of the CRB as shown in papers [40, 41]. It is only in the cases where the MPCs heavily overlap, specially with the LOS, in the directional and time domains that the negative terms are dominant, and then the presence of MPCs degrades the MS localization.

IV Position and Orientation Estimation: Estimator in Beamspace

Next, we propose the use of a beamspace channel transformation in order to estimate the channel parameters in (6). The considered beamspace representation of the channel reduces the complexity by exploiting the sparsity of the mm-wave MIMO channel. If the fractional bandwidth and the number of antennas are not violating the condition for the small array dispersion [11], there exists a common sparse support across all subcarriers. Consequently, the DCS-SOMP method from [8] can be applied for the estimation of AOA, AOD, and TOA. As the estimates of AOA and AOD are limited to lie on a grid defined by the transformation, we apply a refinement of the estimates of all parameters using the SAGE algorithm. Finally, we invoke the EXIP to solve for the position $\mathbf{p}$ and orientation $\alpha$.

IV-A Beamspace Channel Representation

We introduce the $N_{t}\times N_{t}$ transformation matrix, uniformly sampling the virtual spatial angles [44]

[TABLE]

where we assumed $N_{t}$ to be even. Similarly, we define the $N_{r}\times N_{r}$ matrix $\mathbf{U}_{\mathrm{Rx}}$. Both $\mathbf{U}_{\mathrm{Tx}}$ and $\mathbf{U}_{\mathrm{Rx}}$ are unitary matrices. The partial virtual representation of the channel with respect to the angular domain can be written as

[TABLE]

It is readily verified that [11]

[TABLE]

for $-(N_{r}-1)/2\leq i\leq(N_{r}-1)/2$ and $-(N_{t}-1)/2\leq i^{\prime}\leq(N_{t}-1)/2$. We have introduced

[TABLE]

From (33), it is observed that $\check{\mathbf{H}}[n]$ is approximately sparse, since ‘strong’ components are only present in the directions of $\{\theta_{\mathrm{Tx},k}\}$ and $\{\theta_{\mathrm{Rx},k}\}$.

Stacking the observation $\mathbf{y}^{(g)}[n]$ from (6), we obtain

[TABLE]

where

[TABLE]

Hence, since $\check{\mathbf{h}}[n]$ is an approximately sparse vector, we can interpret solving (36) for $\check{\mathbf{h}}[n]$ as a CS problem, allowing us to utilize tools from that domain. In principle, the columns of $\mathbf{U}_{\mathrm{Tx}}$ and $\mathbf{U}_{\mathrm{Rx}}$ corresponding to non-zero entries of the sparse vector $\check{\mathbf{h}}[n]$ correspond to coarse estimates of the AOA/AOD, while the entries in $\check{\mathbf{h}}[n]$ are estimates of $\gamma_{n}(h_{k},\tau_{k})$ (including the effect of the functions $\chi_{t}(\cdot)$ and $\chi_{r}(\cdot)$). The latter values can then be used to estimate $\tau_{k}$ for each path. Since the vectors $\check{\mathbf{h}}[n]\in\mathbb{C}^{N_{r}N_{t}\times 1}$, for $i=1,\ldots,N$, corresponding to the sensing matrix $\mathbf{\Omega}[n]$ in (36) are approximately jointly $(K+1)$-sparse, i.e., the support of $\check{\mathbf{h}}[n]$ does not vary significantly from subcarrier to subcarrier, we can use specialized techniques, such as DCS-SOMP for estimating all $\check{\mathbf{h}}[n]$ jointly in an efficient manner.

Based on the above discussion, we propose to use the following approach:

Coarse estimation of AOA/AOD using a modified DCS-SOMP algorithm. 2. 2.

Fine estimation using the SAGE algorithm, initialized by the coarse estimates. 3. 3.

Estimation of the position and orientation.

Remark

The above sparse representation is not unique. Another representation could rely on a sparse vector of length $N_{t}\times N_{r}\times N$, where each entry would then correspond to an AOA/AOD/TOA triplet. However, the complexity of such an approach would be significantly higher, since $N$ is generally a large number.

IV-B Step 1: Coarse Estimation of Channel Parameters using DCS-SOMP

The first stage of the algorithm involves calling the DCS-SOMP algorithm, providing estimates of the number of paths, the AOA/AOD, and estimates of $\check{\mathbf{h}}[n]$. For the sake of completeness, the steps of DCS-SOMP can be found in Algorithm 1. We note that the algorithm is rank-blind as it does not assume knowledge of the number of the paths (i.e., $K+1$) [45]. Since $K+1$ is unknown, we use the change of residual fitting error $\sum_{n=0}^{N-1}\|\mathbf{r}_{t-1}[n]-\mathbf{r}_{t-2}[n]\|^{2}_{2}$ at each iteration $t$ to a threshold $\delta$. The value for $\delta$ is obtained using a similar procedure as in [15]:

[TABLE]

in which $\gamma^{-1}\left(N,x\right)$ denotes the inverse of the incomplete gamma distribution, $\Gamma(N)$ is the gamma function, and ${P}_{\mathrm{fa}}$ is the false alarm probability.

For each path $k=0,\ldots,\hat{K}$, we can now write

[TABLE]

where $\hat{\check{\mathbf{h}}}^{(k)}=[\hat{\check{h}}^{(k)}[0],\ldots,\hat{\check{h}}^{(k)}[N-1]]^{\mathrm{T}}$ in which $\hat{\check{h}}^{(k)}[n]$ is the entry on subcarrier $n$, related to the $k$-th path found in Algorithm 1, $\mathbf{A}(\tau_{k})=\mathrm{diag}\{1,\ldots,e^{-j2\pi(N-1)\tau_{k}/(NT_{s})}\}$, $\mathbf{v}_{k}$ is the $N\times 1$ noise vector, and $\mathbf{z}^{(k)}$ has entries

[TABLE]

For the purpose of coarse estimation, we ignore the dependence on $n$ in (50), leading to the simple model

[TABLE]

where $\mathbf{a}(\tau_{k})=[1,\ldots,e^{-j2\pi(N-1)\tau_{k}/(NT_{s})}]^{\mathrm{T}}$ and ${z}^{(k)}$ is as in (50), but considering only $\lambda_{c}$ instead of $\lambda_{n}$. From this model, we can recover $\tau_{k}$ and $\tilde{h}_{k}$ by solving a least squares (LS) problem

[TABLE]

Solving for ${\tilde{h}}_{k}$ yields

[TABLE]

Substituting (53) into (52) and expanding the square allows us to solve for ${\tau}_{k}$:

[TABLE]

IV-C Step 2: Fine Estimation of Channel Parameters using SAGE

Channel parameter estimates are refined in an iterative procedure, which is initialized by the estimates from step 1. In principle, we can perform an iterative ascent algorithm directly on the log-likelihood function associated with the model (36). However, this requires a multi-dimensional minimization and computationally complex solutions. A more practical approach is to use the SAGE algorithm with the incomplete data space in (36) as the superposition of $K+1$ complete data space $\check{\mathbf{y}}_{k}[n]$ as:

[TABLE]

where $\check{\mathbf{h}}_{k}[n]$ denotes the vectorized form of $\check{\mathbf{H}}_{k}[n]=\mathbf{U}_{\mathrm{Rx}}^{\mathrm{H}}\mathbf{H}_{k}[n]\mathbf{U}_{\mathrm{Tx}}$ with $\mathbf{H}_{k}[n]$ being the corresponding term for the $k$-th path in the channel frequency response $\mathbf{H}[n]$ in (1). Writing (55) for all the subcarriers results in:

[TABLE]

where

[TABLE]

In the $(m+1)$-th iteration where $m$ is the iteration index, the expectation and maximization steps are performed as described below. For the initialization of the iterative procedure, we use the AOA/AOD, TOA, and channel coefficients from the detection phase using $\hat{\theta}^{(0)}_{\mathrm{Tx},k}$ and $\hat{\theta}^{(0)}_{\mathrm{Rx},k}$ obtained from (44) and (45), respectively, $\hat{\tau}^{(0)}_{k}$ computed from (54), and the corresponding coefficient obtained from (53).

Expectation step

We compute the conditional expectation of the hidden data space $\check{\mathbf{y}}_{k}$ log-likelihood function based on the previous estimation $\hat{\boldsymbol{\eta}}^{(m)}$ and the incomplete data space $\check{\mathbf{y}}$ as:

[TABLE]

For $k=0,\ldots,\hat{K}$, we obtain

[TABLE]

where $\check{\boldsymbol{\mu}}(\boldsymbol{\eta}_{k})=\check{\mathbf{\Omega}}\check{\mathbf{h}}_{k}$, and

[TABLE]

Maximization step

The goal is to find $\boldsymbol{\eta}_{k}$ such that (58) is maximized. In other words, we have

[TABLE]

Solving (60) directly for $\boldsymbol{\eta}_{k}$ is analytically complex due to the fact that it is hard to compute the gradient and Hessian with respect to $\boldsymbol{\eta}_{k}$. Instead, we update the parameters $\hat{\theta}^{(m+1)}_{\mathrm{Tx},k}$, $\hat{\theta}^{(m+1)}_{\mathrm{Rx},k}$, $\hat{\tau}^{(m+1)}_{k}$, and $\hat{\tilde{h}}^{(m+1)}_{k}$ sequentially using Gauss-Seidel-type iterations [46].

IV-D Step 3: Conversion to Position and Rotation Angle Estimates

As a final step, based on the refined estimates of AOA/AOD/TOA from step 2, here we show how the position and orientation of the MS is recovered. Four scenarios are considered: LOS, NLOS, OLOS, and unknown condition.

• •

LOS: When $\hat{K}=1$ and we are in LOS condition, the expressions (17), (19), and (22) describe a mapping $\boldsymbol{\eta}=\boldsymbol{f}_{\mathrm{los}}(\tilde{\boldsymbol{\eta}})$. The classical invariance principle of estimation theory is invoked to prove the equivalence of minimizing the maximum likelihood (ML) criterion in terms of either $\boldsymbol{\eta}_{0}$ or $\tilde{\boldsymbol{\eta}}_{0}$ [47]. Consequently, the estimated values of $\hat{\mathbf{p}}$ and $\hat{\alpha}$ are obtained directly from

[TABLE]

• •

NLOS: For the case with $\hat{K}$ scatterers and a LOS path, the EXIP can be used, as (17)–(22) describe a mapping $\boldsymbol{\eta}=\boldsymbol{f}_{\mathrm{nlos}}(\tilde{\boldsymbol{\eta}})$. Consequently, the estimated $\hat{\tilde{\boldsymbol{\eta}}}$ obtained as

[TABLE]

is asymptotically (w.r.t. $G\times N$) equivalent to the ML estimate of the transformed parameter $\tilde{\boldsymbol{\eta}}$ [35, 36]. Note that $\mathbf{J}_{\boldsymbol{\eta}}$ could be replaced by the identity matrix, leading also to a meaningful estimator of $\tilde{\boldsymbol{\eta}}$, although with probably slightly larger root-mean-square error (RMSE). The Levenberg-Marquardt algorithm (LMA) can be used to solve (63) [48, 49], initialized as follows: we first estimate $\hat{\mathbf{p}}$ and $\hat{\alpha}$ from the LOS path (i.e., the path with the smallest delay). Then, for the first-order reflection $\hat{\mathbf{s}}_{k}$ can be obtained by the intersection of the following two lines: $\tan(\pi-(\hat{\theta}_{\mathrm{Rx},k}+\hat{\alpha}))=(\hat{p}_{y}-s_{1,y})/(\hat{p}_{x}-s_{1,x})$ and $\tan(\hat{\theta}_{\mathrm{Tx},k})=(s_{1,y}-q_{y})/(s_{1,x}-q_{x})$.

• •

OLOS: For the case with $\hat{K}$ scatterers and no LOS path, the EXIP could be used, as (18), (20), and (21) describe a mapping $\boldsymbol{\eta}_{\mathrm{olos}}=\boldsymbol{f}_{\mathrm{olos}}(\tilde{\boldsymbol{\eta}}_{\mathrm{olos}})$. Consequently, the estimated $\hat{\tilde{\boldsymbol{\eta}}}_{\mathrm{olos}}$ obtained as

[TABLE]

is asymptotically equivalent to the ML estimate of the transformed parameter $\tilde{\boldsymbol{\eta}}_{\mathrm{olos}}$ where $\mathbf{J}_{\hat{\boldsymbol{\eta}}_{\mathrm{olos}}}$ denotes the FIM of $\boldsymbol{\eta}_{\mathrm{olos}}$. The estimated parameters from the NLOS links could be used to initialize $\tilde{\boldsymbol{\eta}}_{\mathrm{olos}}$ for the application of the LMA algorithm. The process is slightly more involved than under NLOS. We consider different trial values of $\alpha$, with a resolution $\Delta\alpha$ over a range $[-\alpha_{m},+\alpha_{m}]$ of possible rotation values. For each trial value $\hat{\alpha}_{\mathrm{trial}}$, we can find a corresponding estimate of $\mathbf{p}$. For instance, by solving a set of linear equations for two paths:

[TABLE]

where $d_{k,1}$ was introduced in Fig. 1. After solving (65) for $[\mathbf{p},d_{1,1},d_{2,1}]$, it is straightforward to determine the scatterer locations (as was done in the NLOS case). For each trial value $\hat{\alpha}_{\mathrm{trial}}$, we can then apply the LMA to (64) to obtain $\hat{\tilde{\boldsymbol{\eta}}}_{\mathrm{olos}}$. The solution $\hat{\tilde{\boldsymbol{\eta}}}_{\mathrm{olos}}$ with the smallest $v_{\mathrm{olos}}(\tilde{\boldsymbol{\eta}}_{\mathrm{olos}})$ (with respect to all possible trial value $\hat{\alpha}_{\mathrm{trial}}$) is then retained. Clearly, there is a performance/complexity trade-off based on the choice of $\Delta\alpha$. It is readily seen that to obtain estimates of all parameters, at least three scatterers are needed, since then we have 9 available estimated parameters (1 AOA, 1 AOD, 1 TOA per path) and 9 unknowns (6 scalars for the scatterer locations $\mathbf{s}_{k}$, 3 scalars for $\mathbf{p}$ and $\alpha$).

• •

Unknown: For the case that the receiver does not know whether it operates in NLOS or OLOS, the receiver could apply the technique above under NLOS and under OLOS, separately. This will give two solutions with different cost (measured in terms of (63) and (64)). The best solution (the one with lowest cost) can then be retained.

The complexity analysis for each step of the aforementioned algorithm is presented in Appendix B.

V Simulation Results

In this section, we present simulation results show the values of the bounds and the performance of the proposed estimators for different parameters.

V-A Simulation Setup

We consider a scenario representative of indoor localization in a small conference room with the maximum distance between MS and BS of 4 meters [50]. We set $f_{c}=60\>\mathrm{GHz}$, $B=100\>\mathrm{MHz}$, $c=0.299792\>\mathrm{m}/\mathrm{ns}$, and $N=20$. The geometry-based statistical path loss is used with path length $d_{k}$ and the number of reflectors in each path is set to one, i.e., it is assumed that there is one reflector in each NLOS path [51]. The path loss $\rho_{k}$ between BS and MS for the $k$-th path is computed based on geometry statistics [52, 53]. We set

[TABLE]

where $\sigma^{2}_{0}$ is the reflection loss, $\mathbb{P}_{0}(d_{k,2})=(\gamma_{r}d_{k,2})^{2}e^{-\gamma_{r}d_{k,2}}$ denotes the Poisson distribution of environment geometry with density $\gamma_{r}$ (set to $1/7$ [51]), $\xi^{2}(d_{k})$ denotes the atmospheric attenuation over distance $d_{k}$, and the last term is the free space path loss over distance $d_{k}$. For the LOS link, we obtain

[TABLE]

The average reflection loss for the first-order reflection $\sigma^{2}_{0}$ is set to $-10$ dB with the root-mean-square (RMS) deviation equal to $4$ dB [54], and the atmospheric attenuation over distance $d_{k}$ is set to $16$ dB/km [3]. The number of transmit and receive antennas are set to $N_{t}=65$ and $N_{r}=65$, respectively. The number of simultaneous beams is $M_{t}=1$, and the number of sequentially transmitted signals is $G=32$, unless otherwise stated. The BS is located at $\mathbf{q}\>[\mathrm{m}]=[0,0]^{\mathrm{T}}$ and the MS is located at $\mathbf{p}\>[\mathrm{m}]=[4,0]^{\mathrm{T}}$ with the rotation angle $\alpha=0.1\>\mathrm{rad}$. The elements of the analog beamformers are generated as random values uniformly distributed on the unit circle. The sequences $\tilde{\mathbf{x}}^{(g)}[n]=\mathbf{F}^{(g)}_{\mathrm{BB}}[n]\mathbf{x}^{(g)}[n]$ are obtained as complex exponential terms $e^{j\phi_{g,n}}$ with uniform random phases in $[0,2\pi)$ along different subcarriers, indexed by $n$, and sequentially transmitted symbols, indexed by $g$. The values of the CRB for $\sqrt{\mathrm{CRB}(\tau_{k})}$, $\sqrt{\mathrm{CRB}(\theta_{\mathrm{Rx},k})}$, and $\sqrt{\mathrm{CRB}(\theta_{\mathrm{Tx},k})}$ are defined similar to PEB and REB in (27) and (28), that is, by inverting the FIM $\mathbf{J}_{\tilde{\boldsymbol{\eta}}}$ from (15), choosing the corresponding diagonal entries and taking the square root. Finally, the received signal-to-noise ratio (SNR) is defined as

[TABLE]

in which $\mathrm{diag}\{\cdot\}$ creates a block diagonal matrix from its arguments and $\mathrm{vec}\{\cdot\}$ creates a tall column vector from its arguments.

The performance of the RMSE of the estimation algorithm was assessed from $1000$ Monte Carlo realizations. The false alarm probability was set to ${P}_{\mathrm{fa}}=10^{-3}$ to determine the threshold $\delta$.

V-B Results and Discussion

The Performance versus number of sequential beams

Fig. 2 shows the cumulative distribution function (CDF) of the PEB and the $\mathrm{RMSE}(\hat{\mathbf{p}})$ as a function of the number of beams for LOS conditions. The MS can be anywhere in a rectangle with vertices at the coordinates (in meters): $(2,0)$, $(4,0)$, $(2,0.3)$, and $(4,0.3)$. The signal is scaled so that the total transmit power is kept constant. By increasing the number of beams $G$, the probability of covering the target location in the specified area with a certain accuracy increases. In other words, due to the ergodicity of the process localization accuracy with a certain number of randomly selected sequential beams in each step converges to a constant value for sufficient number of beams $G$. The reason is that for a larger number of beams, the bound decreases thanks to the better spatial coverage. But this effect vanished when the number of beams is sufficient to cover the area where the MS may be located, and then increasing the number of beams only translated into an increased complexity. In principle, the 3 dB beam width for the ULA is approximately $2/N_{t}$, thus reducing when increasing the number of transmit antennas $N_{t}$. Consequently, the number of required beams $G$ to cover the target location in the specified area with the same probability increases. Similarly, by reducing the number of transmit antennas $N_{t}$, the number of required beams $G$ to cover the area decreases. However, the localization accuracy is improved for the case with larger number of transmit antennas $N_{t}$ with the cost of transmitting more beams $G$ for the same coverage. It is observed that for the aforementioned system parameters, $G\geq 20$ randomly selected beams approximately provides the same localization accuracy with $\mathrm{CDF}=0.9$. Note that fewer beams would be needed under a well-chosen deterministic strategy. The same behavior has been observed in NLOS conditions.

Performance in LOS

Fig. 3 shows the evolution of the RMSE of TOA and AOA/AOD in the LOS conditions. The Cramér-Rao bounds are shown by the red lines with the corresponding markers. It is observed that after a few iterations of Algorithm 2 the RMSE of TOA and AOA/AOD converges to the corresponding bounds even for $\mathrm{SNR}=-20\>\mathrm{dB},-10\>\mathrm{dB},0\>\mathrm{dB}$. The performance of the RMSE of the estimation algorithm with respect to different values of the received SNR is shown in Fig. 4–5. It is observed that after $\mathrm{SNR}\approx-20$ dB the RMSE of the TOA, AOA/AOD, rotation angle, and position converge to their corresponding bounds (red dashed lines). Moreover, the proposed algorithm performs well even for very low values of the received SNR, which is the typical case at mm-wave systems before beamforming. We observe that at $\mathrm{SNR}\approx-20\>\mathrm{dB}$ the TOA, AOA/AOD, rotation angle, and position approach the corresponding bounds.

Performance in NLOS

Fig. 6 shows the evolution of the RMSE of TOA and AOA/AOD for $1000$ Monte Carlo realizations in the presence of a scatterer located at $\mathbf{s}_{k}\>[\mathrm{m}]=[1.5,0.4]^{\mathrm{T}}$. It can be observed that the RMSE of the TOA and the AOA/AOD obtained with the proposed algorithm for both the parameters of the LOS and the reflected signals converges to the theoretical also in this case, even at very low received SNR. At $\mathrm{SNR}\approx-5\>\mathrm{dB}$ the TOA, AOA/AOD, rotation angle, and position approach the corresponding bounds.

Performance in OLOS

Finally, the performance in the OLOS case for three scatterers located at $\mathbf{s}_{k}\>[\mathrm{m}]=[1.5,0.4+0.5(k-1)]^{\mathrm{T}}$ for $k=1,2,3$ is investigated in this section using two different initializations of the rotation angle: one with grid resolution $\Delta\alpha\>[\mathrm{rad}]=0.01$ and one with $\Delta\alpha\>[\mathrm{rad}]=0.05$. For both, we set $\alpha_{m}\>[\mathrm{rad}]=0.5$. Fig. 9 shows the performance of the RMSE with respect to the received SNR for position and rotation angle estimation. The proposed estimation method approaches the bound even for the initialization with the resolution $\Delta\alpha\>[\mathrm{rad}]=0.05$. However, the performance of the estimation algorithm is dependent on the resolution of the grid of points $\Delta\alpha$. In particular, a finer grid for the rotation angle leads to better initial estimates and thus a lower final RMSE. For $\mathrm{SNR}\approx-10\>\mathrm{dB}$ the RMSE of position and rotation angle approach the corresponding bounds. We note that the OLOS values, for a fixed SNR, are significantly higher in the OLOS than in the NLOS case.

Unknown Conditions

To analyze the application of the algorithm when the propagation conditions are unknown, we consider the case where there are three scatterers and the LOS path is blocked, that is, the OLOS condition. Starting with the wrong assumption that the path with the shortest delay is the LOS path (i.e., the NLOS condition) leads to very large values of the cost function (63) compared to the actual value of the cost function (64). The results are summarized in Table. I for the average value of the ratio ${\Delta v}\triangleq v_{\mathrm{nlos}}(\hat{\tilde{\boldsymbol{\eta}}})/v_{\mathrm{olos}}(\hat{\tilde{\boldsymbol{\eta}}}_{\mathrm{olos}})$ between the cost function with the wrong and true assumptions. The values in Table I are obtained by averaging 100 realizations, and with a grid resolution of $\Delta\alpha=0.05[\mathrm{rad}]$. The slight difference in the ratio for different values of SNR is due to the limited number of trials.

It is clear that using the wrong assumption about the path with the shortest delay leads to much larger values of the cost function, i.e., the mean value of the ratio ${\Delta v}$ between the cost function with the wrong and true assumptions is on the order of $5$. The main reason for the increase of the cost function using the wrong assumption about the shortest path is that the estimate of MS rotation angle obtained from the AOA and AOD of this path is heavily erroneous. When the shortest path is considered to be a LOS but it is really a reflection, there is a clear mismatch between the geometry of the propagation and the model equations, since there is a scatterer that breaks the direct relation between AOA and AOD existing with the LOS. This mismatch causes a large error in the initial position that is propagated to the final solution. Therefore, observing the ratio of cost functions, we can identify that the path with the shortest delay is related to the scatterer and the LOS path does not exist, that is to say, the OLOS condition is correctly recognized.

Comparison of LOS versus NLOS Performance

Fig. 10 compares the performance of the positioning algorithm in LOS and NLOS for $\mathrm{SNR}=-5$ dB and $G=20$. The MS is anywhere in the same rectangle described at the beginning of Sec. V-B. The scatterers are located at coordinates (in meters) $\mathbf{s}_{1}=(1.5,0.4)$ and $\mathbf{s}_{2}=(1.5,0.6)$. The accuracy and robustness of the localization algorithm is improved by adding the scatterers compared to the case when only LOS is used. Moreover, the performance in the OLOS is much worse than in LOS or NLOS due to the severe effect of path loss as shown already in the paper by comparing Figs. 5 and 8 with Fig. 9.

VI Conclusion

We have studied the determination of a receiver position and orientation using a single transmitter in a MIMO system. Our study includes LOS, as well as NLOS and OLOS conditions, shedding insight into the potential of locating a receiver even when the LOS is blocked. We have derived fundamental performance bounds on the estimation uncertainty for delay, angle of arrival, angle of departure, and channel gain of each path, as well as the user position and orientation angle. We also proposed a novel three stage algorithm for the estimation of the user position and orientation angle. This algorithm determines coarse estimates of the channel parameters by exploiting the sparsity of the mm-wave in beamspace, followed by an iterative refinement, and finally a conversion to position and orientation. Through simulation studies, we demonstrate the efficiency of the proposed algorithm, and show that even in OLOS conditions, it is possible to estimate the user’s position and orientation angle, by exploiting the information coming from the multipath, though at a significant performance penalty.

Appendix A Elements in (14)

Replacing $\mathbf{y}[n]$ from (6) in (11), using (13), and considering $\mathbb{E}_{\mathbf{y}|\boldsymbol{\eta}}[\mathbf{n}[n]]=\mathbf{0}$, we obtain

[TABLE]

The elements of the FIM are obtained based on (69). The entry associated with the $n$-th subcarrier is denoted as $\Psi_{n}(x_{r},x_{s})$, and given by (for $\left\{\tau_{r},\tau_{s}\right\}$ and $\left\{\boldsymbol{\theta}_{r},\boldsymbol{\theta}_{s}\right\}$)

[TABLE]

The following notations are introduced:

[TABLE]

where $l\in\{0,1\}$, and $\mathbf{A}_{k,n}(\tau_{r},\tau_{s}),k\in\{0,1,2\}$, is given by

[TABLE]

The vectors $\mathbf{a}_{\mathrm{Tx},\mathbf{F},n}(\theta_{\mathrm{Tx},r})$ and $\mathbf{a}_{\mathbf{D}_{\mathrm{Tx}},\mathbf{F},n}(\theta_{\mathrm{Tx},r})$ are given by $\mathbf{a}_{\mathrm{Tx},\mathbf{F},n}(\theta_{\mathrm{Tx},r})=\mathbf{F}^{\mathrm{H}}[n]\mathbf{a}_{\mathrm{Tx},n}(\theta_{\mathrm{Tx},r})$ and $\mathbf{a}_{\mathbf{D}_{\mathrm{Tx}},\mathbf{F},n}(\theta_{\mathrm{Tx},r})=\mathbf{F}^{\mathrm{H}}[n]\mathbf{D}_{\mathrm{Tx},r}[n]\mathbf{a}_{\mathrm{Tx},n}(\theta_{\mathrm{Tx},r})$. The matrix $\mathbf{D}_{\mathrm{Tx},r}[n]$ is defined as

[TABLE]

The scalars $A_{\mathrm{Rx},n}(\theta_{\mathrm{Rx,r}},\theta_{\mathrm{Rx,s}})$, $A_{\mathbf{D}_{\mathrm{Rx,s}},n}(\theta_{\mathrm{Rx,r}},\theta_{\mathrm{Rx,s}})$, and $A_{\mathbf{D}_{\mathrm{Rx,r,s}},n}(\theta_{\mathrm{Rx,r}},\theta_{\mathrm{Rx,s}})$ are defined as

[TABLE]

where $\mathbf{D}_{\mathrm{Rx},r}[n]$ has the same expression as (81) by replacing the subscript $\mathrm{Tx}$ by $\mathrm{Rx}$ and $N_{t}$ by $N_{r}$. The terms including channel coefficients are summarized as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Complexity Analysis

We analyze the complexity of different stages of the proposed algorithm.

• •

Coarse Estimation: The complexity in performing (42) is on the order of $O(N^{2}_{r}N^{2}_{t}GN_{\mathrm{sub}})$ where $N_{\mathrm{sub}}$ denotes the few subcarriers sufficient to detect the dominant path. The QR factorization of the mutilated basis $\mathbf{\Omega}_{\mathcal{K}_{t}}[n]$ approximately requires $O(GN_{r}\hat{K}^{2})$ operations for each subcarrier, and matrix inversion to obtain the channel coefficients in (48) approximately takes $O(N\hat{K}^{3})$ operations for all the subcarriers. The complexity in computing (54) is on the order of $O(ND_{o}\hat{K})$ where $D_{o}$ denotes the number of delay grid points, and (53) requires $O(N\hat{K})$ operations. Consequently, the maximum complexity from coarse estimation of the channel parameters is dominated by the term $\hat{K}\times O(N^{2}_{r}N^{2}_{t}GN_{\mathrm{sub}})$.

• •

Fine Estimation: In the refinement phase, the complexity is mainly affected by Gauss-Seidel-type iterations with first and second order derivatives of a vector $\mathbf{a}(x)$ of length $L_{\mathrm{x}}$ with respect to a variable $x$ that can be delay, AOA, and AOD. These operations lead to a complexity on the order of $O(L^{2}_{\mathrm{x}}N)$ for each path. Given the subsequent path refinement, the maximum complexity of fine estimation is on the order of $O(\hat{K}^{2})\times O(L^{2}_{\mathrm{x}}N)$.

• •

Conversion to Position and Orientation: The conversion to position and orientation in the LOS case is easy to implement since it involves only some basic operations. For the NLOS and OLOS scenarios, the LMA algorithm is applied. It is not considered the complexity driver, since it combines the advantages of gradient-descent and Gauss-Newton methods. The LMA algorithm can be effectively applied by implementing delayed gratification, which leads to higher success rate and fewer Jacobian evaluations.