$\alpha(M,M+\epsilon,M+\epsilon,\mathsf{r}^{\prime},\mathsf{s}^{\prime})+\beta(\mathsf{r}^{\prime},\mathsf{L}^{\prime})\leq\delta(\epsilon).$

This ensures that the auxiliary system (64) satisfies
$\|(\mathsf{r}\mathsf{s}^{2}/\mathsf{L})\cdot\overline{u}(t)-u(t)\|\leq\epsilon$,
provided that Assumption 1 holds. Then, it follows that (55) holds, so that
analogously to the proofs of Proposition 4 and Theorem 1, the condition (56)
guarantees that (60) holds for all $t\geq 0$. It follows that

$g(\mathsf{Dec}(\mathbf{u}(t)))=g(\tilde{\mathsf{z}}(t))=\frac{\mathsf{r}\mathsf{s}^{2}}{\mathsf{L}}\cdot\overline{u}(t)$

holds for all $t\geq 0$, so the proof is completed. $\blacksquare$

Finally, we review the implications of Theorem 2. There have been two issues
that hinder the unlimited operation of linear dynamic systems implemented over
encrypted data; the first issue was due to the recursive multiplication by
non-integer numbers, and it has been handled by the introduced method of
converting the state matrix to integers.

Then, the second issue was the “error growth” problem of the cryptosystems.
Even though the given state matrix of the system consists of only integers,
most homomorphic cryptosystems (that allow for both the addition and
multiplication) do not support the recursive multiplication by the encrypted
numbers for an infinite number of times, unless bootstrapping techniques of
fully homomorphic encryption is utilized. And, this is because of the error
growth problem; once the injected errors exceed a certain bound, then the
correct computation outcome cannot be expected.

To allow for the unlimited recursive multiplication, the benefit of GSW-LWE
cryptosystem has been first discussed. Considering that linear systems only
update the state (multiplicand) and does not update the matrix (multiplier),
it has been observed that the matrix encrypted using GSW scheme can be
multiplied to the LWE-based encrypted vectors for an infinite number of times,
thanks to the structure of LWE-type ciphertexts where the errors are stored
together with the messages. And, despite the recursive multiplication repeated
unlimited times, it has been shown that the growth of the injected errors does
not tend to be infinitely large, and its effect is suppressed and controlled,
under stability. The effect of the errors can be identified with perturbations
(or external disturbances), and its size can be made small by appropriate
choice of parameters, so that it can keep the performance as practically
equivalent to the given un-encrypted model.

Compared to the methods from cryptography utilizing the bootstrapping
techniques of fully homomorphic encryptions, the introduced method exploits
addition and multiplication over ciphertexts only, which does not require a
substantial amount of computational resources. A couple of related results are
introduced; a guideline of choosing the cryptosystem parameters
$(n,q,\sigma_{0})$ and the control parameters
$(\mathsf{r},\mathsf{s},\mathsf{L})$ that guarantees both the desired level of
security and control performance is found in [62]. And, further illustrative
explanations using MATLAB example codes can be found in [63].

##  4 Conclusion

We have introduced several approaches for encrypted control, based on
homomorphic encryption, multi-party computation, and secret sharing. We have
compared their benefits and weaknesses, and discussed the trade-offs between
them. It has been suggested that the encrypted control approach and the
corresponding cryptosystem should be chosen with an engineering point of view,
so that its security model, enabled operations, and computation efficiency is
suitable for the system under construction.

Among the research directions struggling to overcome the trade-offs between
security, computation efficiency, and the operation ability, we have
introduced a homomorphic encryption based problem which aims for implementing
linear dynamic systems exploiting only addition and multiplication over
integers. Then, it has been addressed that the problems from cryptography can
be tackled with the control perspectives; we have addressed that the error
growth of the GSW-LWE cryptosystem can be regarded as perturbations in control
systems, and the issue of incapability of recursive multiplication by non-
integer state matrix has been resolved using a method based on pole-placement
technique.

##  5 Acknowledgement

This study was supported in part by the Research Program funded by the
SeoulTech(Seoul National University of Science and Technology), and in part by
the Knut and Alice Wallenberg Foundation.

## References

  * [1] H. Sandberg, S. Amin, K. H. Johansson, Cyberphysical security in networked control systems: An introduction to the issue, IEEE Control Systems Magazine 35 (1) (2015) 20–23.
  * [2] S. Amin, A. A. Cárdenas, S. S. Sastry, Safe and secure networked control systems under denial-of-service attacks, in: International Workshop on Hybrid Systems: Computation and Control, Springer, 2009, pp. 31–45.
  * [3] R. Langner, Stuxnet: Dissecting a cyberwarfare weapon, IEEE Security & Privacy 9 (3) (2011) 49–51.
  * [4] Y. Liu, P. Ning, M. K. Reiter, False data injection attacks against state estimation in electric power grids, ACM Transactions on Information and System Security (TISSEC) 14 (1) (2011) 1–33.
  * [5] J. Slay, M. Miller, Lessons learned from the maroochy water breach, in: International conference on critical infrastructure protection, Springer, 2007, pp. 73–82.
  * [6] A. Teixeira, I. Shames, H. Sandberg, K. H. Johansson, A secure control framework for resource-limited adversaries, Automatica 51 (2015) 135–148.
  * [7] F. Pasqualetti, F. Dörfler, F. Bullo, Attack detection and identification in cyber-physical systems, IEEE transactions on automatic control 58 (11) (2013) 2715–2729.
  * [8] Y. Mo, B. Sinopoli, Secure control against replay attacks, in: 2009 47th annual Allerton conference on communication, control, and computing (Allerton), IEEE, 2009, pp. 911–918.
  * [9] R. S. Smith, Covert misappropriation of networked control systems: Presenting a feedback structure, IEEE Control Systems Magazine 35 (1) (2015) 82–92.
  * [10] M. Schulze Darup, A. B. Alexandru, D. E. Quevedo, G. J. Pappas, Encrypted control for networked systems: An illustrative introduction and current challenges, IEEE Control Systems Magazine 41 (3) (2021) 58–78.
  * [11] K. Kogiso, T. Fujita, Cyber-security enhancement of networked control systems using homomorphic encryption, in: 2015 54th IEEE Conference on Decision and Control (CDC), IEEE, 2015, pp. 6836–6843.
  * [12] J. Kim, C. Lee, H. Shim, J. H. Cheon, A. Kim, M. Kim, Y. Song, Encrypting controller using fully homomorphic encryption for security of cyber-physical systems, IFAC-PapersOnLine 49 (22) (2016) 175–180.
  * [13] F. Farokhi, I. Shames, N. Batterham, Secure and private control using semi-homomorphic encryption, Control Engineering Practice 67 (2017) 13–20.
  * [14] R. L. Rivest, L. Adleman, M. L. Dertouzos, et al., On data banks and privacy homomorphisms, Foundations of secure computation 4 (11) (1978) 169–180.
  * [15] T. ElGamal, A public key cryptosystem and a signature scheme based on discrete logarithms, IEEE transactions on information theory 31 (4) (1985) 469–472.
  * [16] P. Paillier, Public-key cryptosystems based on composite degree residuosity classes, in: International conference on the theory and applications of cryptographic techniques, Springer, 1999, pp. 223–238.
  * [17] J. H. Cheon, A. Kim, M. Kim, Y. Song, Homomorphic encryption for arithmetic of approximate numbers, in: International Conference on the Theory and Application of Cryptology and Information Security, Springer, 2017, pp. 409–437.
  * [18] C. Dwork, Differential privacy: A survey of results, in: International conference on theory and applications of models of computation, Springer, 2008, pp. 1–19.
  * [19] C. Dwork, A. Roth, et al., The algorithmic foundations of differential privacy., Found. Trends Theor. Comput. Sci. 9 (3-4) (2014) 211–407.
  * [20] O. Regev, On lattices, learning with errors, random linear codes, and cryptography, Journal of the ACM (JACM) 56 (6) (2009) 1–40.
  * [21] V. Lyubashevsky, C. Peikert, O. Regev, On ideal lattices and learning with errors over rings, in: Annual international conference on the theory and applications of cryptographic techniques, Springer, 2010, pp. 1–23.
  * [22] R. Matthews, On the derivation of a “chaotic” encryption algorithm, Cryptologia 13 (1) (1989) 29–42.
  * [23] N. Bourbakis, C. Alexopoulos, Picture data encryption using scan patterns, Pattern Recognition 25 (6) (1992) 567–581.
  * [24] D. Boneh, E.-J. Goh, K. Nissim, Evaluating 2-dnf formulas on ciphertexts, in: Theory of cryptography conference, Springer, 2005, pp. 325–341.
  * [25] M. v. Dijk, C. Gentry, S. Halevi, V. Vaikuntanathan, Fully homomorphic encryption over the integers, in: Annual international conference on the theory and applications of cryptographic techniques, Springer, 2010, pp. 24–43.
  * [26] Z. Brakerski, V. Vaikuntanathan, Fully homomorphic encryption from ring-lwe and security for key dependent messages, in: Annual cryptology conference, Springer, 2011, pp. 505–524.
  * [27] J. Fan, F. Vercauteren, Somewhat practical fully homomorphic encryption, Cryptology ePrint Archive.
  * [28] Z. Brakerski, C. Gentry, V. Vaikuntanathan, (leveled) fully homomorphic encryption without bootstrapping, ACM Transactions on Computation Theory (TOCT) 6 (3) (2014) 1–36.
  * [29] M. Schulze Darup, A. Redder, I. Shames, F. Farokhi, D. Quevedo, Towards encrypted mpc for linear constrained systems, IEEE Control Systems Letters 2 (2) (2017) 195–200.
  * [30] A. B. Alexandru, A. Tsiamis, G. J. Pappas, Towards private data-driven control, in: 2020 59th IEEE Conference on Decision and Control (CDC), IEEE, 2020, pp. 5449–5456.
  * [31] J. Suh, T. Tanaka, Encrypted value iteration and temporal difference learning over leveled homomorphic encryption, in: 2021 American Control Conference (ACC), IEEE, 2021, pp. 2555–2561.
  * [32] C. Gentry, Fully homomorphic encryption using ideal lattices, in: Proceedings of the forty-first annual ACM symposium on Theory of computing, 2009, pp. 169–178.
  * [33] L. Ducas, D. Micciancio, Fhew: bootstrapping homomorphic encryption in less than a second, in: Annual international conference on the theory and applications of cryptographic techniques, Springer, 2015, pp. 617–640.
  * [34] I. Chillotti, N. Gama, M. Georgieva, M. Izabachene, Faster fully homomorphic encryption: Bootstrapping in less than 0.1 seconds, in: international conference on the theory and application of cryptology and information security, Springer, 2016, pp. 3–33.
  * [35] C. Gentry, S. Halevi, N. P. Smart, Better bootstrapping in fully homomorphic encryption, in: International Workshop on Public Key Cryptography, Springer, 2012, pp. 1–16.
  * [36] R. Fritz, M. Fauser, P. Zhang, Controller encryption for discrete event systems, in: 2019 American Control Conference (ACC), IEEE, 2019, pp. 5633–5638.
  * [37] S. Zhou, Z. Yu, E. S. A. Nasr, H. A. Mahmoud, E. M. Awwad, N. Wu, Homomorphic encryption of supervisory control systems using automata, IEEE Access 8 (2020) 147185–147198.
  * [38] A. Shamir, How to share a secret, Communications of the ACM 22 (11) (1979) 612–613.
  * [39] E. Karnin, J. Greene, M. Hellman, On secret sharing systems, IEEE Transactions on Information Theory 29 (1) (1983) 35–41.
  * [40] E. F. Brickell, Some ideal secret sharing schemes, in: Workshop on the Theory and Application of of Cryptographic Techniques, Springer, 1989, pp. 468–475.
  * [41] A. Beimel, Secret-sharing schemes: A survey, in: International conference on coding and cryptology, Springer, 2011, pp. 11–46.
  * [42] M. Schulze Darup, T. Jager, Encrypted cloud-based control using secret sharing with one-time pads, in: 2019 IEEE 58th Conference on Decision and Control (CDC), IEEE, 2019, pp. 7215–7221.
  * [43] Y. Mo, R. M. Murray, Privacy preserving average consensus, IEEE Transactions on Automatic Control 62 (2) (2016) 753–765.
  * [44] E. Nozari, P. Tallapragada, J. Cortés, Differentially private average consensus: Obstructions, trade-offs, and optimal algorithm design, Automatica 81 (2017) 221–231.
  * [45] Y. Wang, Privacy-preserving average consensus via state decomposition, IEEE Transactions on Automatic Control 64 (11) (2019) 4711–4716.
  * [46] S. Schlor, M. Hertneck, S. Wildhagen, F. Allgöwer, Multi-party computation enables secure polynomial control based solely on secret-sharing, in: 2021 60th IEEE Conference on Decision and Control (CDC), IEEE, 2021, pp. 4882–4887.
  * [47] Y. Lindell, How to simulate it–a tutorial on the simulation proof technique, Tutorials on the Foundations of Cryptography (2017) 277–346.
  * [48] I. Damgård, M. Fitzi, E. Kiltz, J. B. Nielsen, T. Toft, Unconditionally secure constant-rounds multi-party computation for equality, comparison, bits and exponentiation, in: Theory of Cryptography Conference, Springer, 2006, pp. 285–304.
  * [49] T. Nishide, K. Ohta, Multiparty computation for interval, equality, and comparison without bit-decomposition protocol, in: International Workshop on Public Key Cryptography, Springer, 2007, pp. 343–360.
  * [50] P. Bogetoft, D. L. Christensen, I. Damgård, M. Geisler, T. Jakobsen, M. Krøigaard, J. D. Nielsen, J. B. Nielsen, K. Nielsen, J. Pagter, et al., Secure multiparty computation goes live, in: International Conference on Financial Cryptography and Data Security, Springer, 2009, pp. 325–343.
  * [51] M. Dahl, C. Ning, T. Toft, On secure two-party integer division, in: International Conference on Financial Cryptography and Data Security, Springer, 2012, pp. 164–178.
  * [52] R. Cramer, I. B. Damgård, et al., Secure multiparty computation, Cambridge University Press, 2015.
  * [53] A. B. Alexandru, M. Morari, G. J. Pappas, Cloud-based mpc with encrypted data, in: 2018 IEEE Conference on Decision and Control (CDC), IEEE, 2018, pp. 5014–5019.
  * [54] A. B. Alexandru, G. J. Pappas, Secure multi-party computation for cloud-based control, in: Privacy in Dynamical Systems, Springer, 2020, pp. 179–207.
  * [55] F. J. Gonzalez-Serrano, A. Amor-Martın, J. Casamayon-Anton, State estimation using an extended kalman filter with privacy-protected observed inputs, in: 2014 IEEE International Workshop on Information Forensics and Security (WIFS), IEEE, 2014, pp. 54–59.
  * [56] C. Gentry, A. Sahai, B. Waters, Homomorphic encryption from learning with errors: Conceptually-simpler, asymptotically-faster, attribute-based, in: Annual Cryptology Conference, Springer, 2013, pp. 75–92.
  * [57] R. Lindner, C. Peikert, Better key sizes (and attacks) for lwe-based encryption, in: Cryptographers’ Track at the RSA Conference, Springer, 2011, pp. 319–339.
  * [58] L. Chen, L. Chen, S. Jordan, Y.-K. Liu, D. Moody, R. Peralta, R. Perlner, D. Smith-Tone, Report on post-quantum cryptography, Vol. 12, US Department of Commerce, National Institute of Standards and Technology, 2016.
  * [59] N. Genise, C. Gentry, S. Halevi, B. Li, D. Micciancio, Homomorphic encryption for finite automata, in: International Conference on the Theory and Application of Cryptology and Information Security, Springer, 2019, pp. 473–502.
  * [60] J. Kim, H. Shim, K. Han, Dynamic controller that operates over homomorphically encrypted data for infinite time horizon, IEEE Transactions on Automatic Control.
  * [61] J. Kim, H. Shim, H. Sandberg, K. H. Johansson, Method for running dynamic systems over encrypted data for infinite time horizon without bootstrapping and re-encryption, in: 2021 60th IEEE Conference on Decision and Control (CDC), IEEE, 2021, pp. 5614–5619.
  * [62] J. Kim, H. Shim, K. Han, Design procedure for dynamic controllers based on lwe-based homomorphic encryption to operate for infinite time horizon, in: 2020 59th IEEE Conference on Decision and Control (CDC), IEEE, 2020, pp. 5463–5468.
  * [63] J. Kim, H. Shim, K. Han, Comprehensive introduction to fully homomorphic encryption for dynamic feedback controller via lwe-based cryptosystem, in: Privacy in Dynamical Systems, Springer, 2020, pp. 209–230.


# Low complexity joint position and channel estimation at millimeter wave
based on multidimensional orthogonal matching pursuit

Joan Palacios, Nuria González-Prelcic
North Carolina State University, USA
<EMAIL_ADDRESS>Cristian Rusu
University of Bucharest, Romania
<EMAIL_ADDRESS>
###### Abstract

Compressive approaches provide a means of effective channel high resolution
channel estimates in millimeter wave MIMO systems, despite the use of analog
and hybrid architectures. Such estimates can also be used as part of a joint
channel estimation and localization solution. Achieving good localization
performance, though, requires high resolution channel estimates and better
methods to exploit those channels. In this paper, we propose a low complexity
multidimensional orthogonal matching pursuit strategy for compressive channel
estimation based by operating with a product of independent dictionaries for
the angular and delay domains, instead of a global large dictionary. This
leads to higher quality channel estimates but with lower complexity than
generalizations of conventional solutions. We couple this new algorithm with a
novel localization formulation that does not rely on the absolute time of
arrival of the LoS path and exploits the structure of reflections in indoor
channels. We show how the new approach is able to operate in realistic 3D
scenarios to estimate the communication channel and locate devices in an
indoor simulation setting.

##  I Introduction

Millimeter wave (mmWave) communication and MIMO technology offer additional
benefits beyond high data rate communications. The large arrays at high
frequencies provide the angle and delay resolvability that enables accurate
localization of users and objects in the environment as a byproduct of
communication [1]. Future 5G releases and 6G envision VR/AR and automated cars
as main use cases, expecting an indoor accuracy $<1$ cm and an outdoor
accuracy $<10$ cm. New research work is needed, from an algorithmic
perspective, to achieve this performance.

One direction for developing high accuracy localization in a mmWave MIMO
system is to take a model-driven two-stage approach [1, 2, 3, 4]. The first
stage performs compressive channel estimation, which has proven successful in
mmWave and massive MIMO channels explicitly, since the large bandwidth or
number of antennas prevents the use of conventional channel estimation
approaches [5, 6, 7, 8]. The second stage exploits the geometric relationships
between some of the parameters of a sparse MIMO channel (angles and delays)
and the position and orientation of the transmitter and the receiver to
complete the localization. Despite the success of compressive approaches,
further work is needed to overcome obstacles to achieving higher localization
accuracy along with relaxing key assumptions that make the algorithms more
practical. To improve accuracy, more precise channel estimators are needed
that do not suffer from high complexity encountered with larger dictionaries.
To improve realizability, synchronization assumptions need to be relaxed. For
example, prior work assumed that the absolute delay of the line of sight (LoS)
component can be directly obtained by the receiver as part of the channel
estimate, which requires an unrealistic system where the transmitter and
receiver are triggered at the same time. To take advantage of larger arrays,
more work is needed on hybrid architectures versus the fully digital
architectures with high resolution converters, which is not feasible at mmWave
frequencies due to its high power consumption [9]. Finally, bandlimited models
are needed to avoid an equivalent filtering effect based on a Dirac delta
function, which leads to an artificial enhancement of the channel sparsity.

In this paper, we present a practical approach to joint channel estimation and
localization at mmWave that overcomes previous limitations. To drastically
reduce complexity, we propose the use of a multidimensional orthogonal
matching pursuit algorithm (MOMP) [10] that operates with a dictionary in
multiple dimensions instead of a large dictionary as conventional OMP [11]
would do. To enable localization with a realistic transceiver, the position
estimate is obtained as a function of the time difference of arrival and the
angular parameters of the different multipath components, without relaying in
an absolute time of arrival of the LoS path provided by the channel estimation
algorithm. In addition, we also propose a novel approach to map channel
parameters into position information by exploiting the special structure of
reflections on walls, ceiling, and floor that appears in indoor channels.
Numerical results in an indoor setting simulated by ray tracing show the
effectiveness of the proposed strategy when the transmitter and receiver
employ a practical hybrid MIMO architecture.

Notation: We use the following notation throughout the paper. $x$, ${\bf x}$,
${\bf X}$ and $\mathcal{X}$ will be the styles for scalar, vector, matrix or
tensor and set. Regarding sub/supper-indices, $x$ and ${\rm x}$ are used to
denote scalar and categorical values respectively. $[{\bf x}]_{n}$ denotes the
$n$-th entry of ${\bf x}$. For a 2D matrix ${\bf X}$, $[{\bf X}]_{a,b}$,
$[{\bf X}]_{a,:}$ and $[{\bf X}]_{:,b}$ are respectively, the element in the
$a$-th row and $b$-th column, the $a$-th row and the $b$-th column, this
notation is extended to the case of tensors with multi-indices acting like
multiple indices, such as $[{\bf X}]_{{\bf a},b}=[{\bf X}]_{a_{1},a_{2},b}$
for ${\bf a}=[a_{1},a_{2}]$. The set $\mathcal{U}^{n\times m}$ denotes the set
of unit magnitude complex matrices of size $n$ by $m$. We use the operator
$\|{\bf x}\|$, $\|{\bf X}\|$ to denote the Euclidean and Frobenius norms of
${\bf x}$ and ${\bf X}$ respectively. In this paper, indexing of tensors
starts from 1 and we consider all vectors are column vectors.

##  II System and signal model

We consider a MIMO communication system with a hybrid architecture at both
ends to enable operation at millimeter wave bands . The number of antennas at
the transmit array is $N_{\rm T}$, while $N_{\rm R}$ denotes the number of
antenna elements at the receive array. The number of RF-chains is denoted as
$M_{\rm T}$ and $M_{\rm R}$ for the transmitter and receiver respectively. We
consider the transmission of $N_{\rm s}$ streams, with $N_{\rm s}\leq M_{\rm
T}$ in general. During the link establishment phase, the transmitter sends a
sequence of $M$ training symbol vectors to the receiver, so this can estimate
the communication channel and the location of the transmitter. The training
symbols are transmitted by a sequence of training hybrid precoders designed to
sound the channel, and received through a set of training hybrid combiners.
For training purposes, we choose $N_{\rm s}=M_{\rm T}$, what leads to a square
digital precoder. The digital combiner is also chosen to be square during
training, to avoid compressing the information captured by the receive RF
chains. With these assumptions in mind, the digital counterpart of the
training precoders/combiners for the $m$-th training frame are ${\bf
F}_{m}^{\rm BB}\in\mathbb{C}^{M_{\rm T}\times M_{\rm T}}$ and ${\bf
W}_{m}^{\rm BB}\in\mathbb{C}^{M_{\rm R}\times M_{\rm R}}$. The analog training
precoder and the analog training combiner for the $m$-th training frame are
denoted as ${\bf F}_{m}^{\rm RF}\in\mathcal{U}^{N_{\rm T}\times M_{\rm T}}$
and ${\bf W}_{m}^{\rm RF}\in\mathcal{U}^{N_{\rm R}\times M_{\rm R}}$. This
way, the hybrid precoders/combiners are ${\bf F}_{m}={\bf F}_{m}^{\rm RF}{\bf
F}_{m}^{\rm BB}\in\mathbb{C}^{N_{\rm T}\times M_{\rm T}}$ and ${\bf
W}_{m}={\bf W}_{m}^{\rm RF}{\bf W}_{m}^{\rm BB}\in\mathbb{C}^{N_{\rm R}\times
M_{\rm R}}$. Regarding the construction of the training sequence itself, we
consider the transmission of $M_{\rm T}$ streams during training, with a
length $D$ zero padding and $Q$ symbols per stream and frame. $D$ is taken as
the delay tap length of the channel. Under these assumptions and definitions,
the training symbol matrix for the $m$-th frame is denoted as
${\mathbf{S}}_{m}\in\mathbb{C}^{M_{\rm T}\times(Q+D)}$.

The frequency selective mmWave channel is modeled using a geometric channel
model with $L$ paths. The $d$-th delay tap of the channel, for $d=1,\ldots,D$,
is represented as

${\bf H}_{d}=\sum_{l=1}^{L}\alpha_{l}{\bf a}_{\rm R}(\bm{\theta}_{l}){\bf
a}_{\rm T}^{\rm H}(\bm{\phi}_{l})p((d-1)T_{\rm s}+\tau_{0}-\tau_{l}),$ (1)

where $\alpha_{l}\in\mathbb{C}$, $\tau_{l}\in\mathbb{R}$, $\bm{\theta}_{l}$,
and $\bm{\phi}_{l}$ are the complex gain, delay, angle of arrival (AoA), and
angle of departure (AoD) for the $l$-th path, $p(t)$ is a band limited pulse
shaping function including filtering effects at the transmitter and receiver,
$\tau_{0}$ is the delay between the beginning of the transmission and the
beginning of the reception (clock offset), and ${\bf a}_{\rm
T}(\bm{\phi}_{l})$ and ${\bf a}_{\rm R}(\bm{\theta}_{l})$ denote the array
response vectors for the transmitter and the receiver. The angular directions
are represented as unitary vector directions, i.e.
${\bm{\theta}}_{l},{\bm{\phi}}_{l}\in\\{{\bf v}\in\mathbb{R}^{3}\text{ such
that }\|{\bf v}\|=1\\}$.

Because of the multiple array configurations exploited simultaneously with a
hybrid architecture, the transmission of each training frame will generate a
set of receive signals that we denote as the block matrix
${{\mathbf{Y}}}_{m}\in\mathbb{C}^{M_{\rm R}\times Q}$, comprised of $M_{\rm
R}$ combinations of the $M_{\rm T}$ pilot signal streams. Considering a
transmission power $P$, the expression of this received block matrix is

$[{{\mathbf{Y}}}_{m}]_{:,q}\\!=\\!\\!\sqrt{P}\sum_{d=1}^{D}\\!{\bf W}_{m}^{\rm
H}{\bf H}_{d}{\bf F}_{m}[{{\mathbf{S}}_{m}}]_{:,q+D-d}\\!+\\!{\bf W}_{m}^{\rm
H}[{{\mathbf{N}}}_{m}]_{:,q},$ (2)

The noise block matrix ${\bf N}_{m}\in\mathbb{C}^{N_{\rm R}\times Q}$ has
independent identically distributed entries following a distribution
$\mathcal{NC}(0,\sigma^{2})$, being $\sigma^{2}$ the noise power. The final
observation matrix ${\mathbf{Y}}$ is constructed by concatenating the received
signals ${{\mathbf{Y}}}_{m}$ for the all the training frames, i.e.
${\mathbf{Y}}=[{\mathbf{Y}}_{1},{\mathbf{Y}}_{2},\ldots,{\mathbf{Y}}_{M}]\in\mathbb{C}^{M_{\rm
R}\times QM}$. The problem to be solved consists of the estimation of the
channel matrices ${\bf H}_{d}$ and the position of the transmitter assuming
that the receiver position is fixed and known.

##  III MOMP-based mmWave channel estimation

Prior work has developed a significant number of sparse recovery solutions to
solve the problem of estimating the mmWave channel . The proposed strategies
incur, however, in high computational complexity. This is mainly due to
exploiting a formulation, for both the sensing matrix and the sparsifying
dictionary, based on Kronecker products of large matrices. The Kronecker
construction leads to an equivalent measurement matrix of significant
dimension for common sizes of mmWave arrays, what makes matching pursuit
solutions impractical in real-world 3D scenarios. In this section, we derive a
new formulation based on MOMP [10], which leverages independent dictionaries
in the angular and delay domains, without building a global dictionary of
larger size.

###  III-A Background in MOMP

The multidimensional matching pursuit problem [10] consists of reconstructing
a multidimensional sparse signal to best fit the available observations,
assuming that it can be represented by projections on a given set of
sparsifying dictionaries. Let us first consider a whitened version of the
observation matrix ${\bf Y}\in\mathbb{C}^{M_{\rm R}\times QM}$ as previously
defined, which contains a set of $QM$ observations, each one of dimension
$M_{\rm R}$. We denote the vectorized whitened measurements as ${\bf
y}\in\mathbb{C}^{MM_{\rm R}Q}$.To reconstruct a sparse tensor from this
observation using MOMP, we consider $N_{\rm D}$ dictionaries, with the $k$-th
dictionary ${\bf\Psi}_{k}\in\mathbb{C}^{N_{k}^{\rm s}\times N_{k}^{\rm a}}$
consisting of $N_{k}^{\rm a}$ atoms in $\mathbb{C}^{N_{k}^{\rm s}}$. The
coefficients of the sparse signal in the set of dictionaries are represented
by ${\bf C}\in\mathbb{C}^{N_{1}^{\rm a}\times\ldots\times N_{N_{\rm D}}^{\rm
a}}$. Our goal is to reconstruct the tensor containing these coefficients
under the assumption that only a few elements are non-zero. To cycle through
the multiple indices of the dictionaries, we define the set of entry
coordinate combinations $\mathcal{I}=\\{{\bf i}=(i_{1},\ldots,i_{N_{\rm
D}})\in\mathbb{N}^{N_{\rm D}}\text{ s.t. }i_{k}\leq N_{k}^{\rm s}\quad\forall
d\leq N_{\rm D}\\}$, and the set of dictionary index combinations
$\mathcal{J}=\\{{\bf j}=(j_{1},\ldots,j_{N_{\rm D}})\in\mathbb{N}^{N_{\rm
D}}\text{ s.t. }j_{k}\leq N_{k}^{\rm a},\forall d\leq N_{\rm D}\\}$. Finally,
the sensing tensor ${\bf\Phi}\in\mathbb{C}^{MM_{\rm R}Q\times N_{1}^{\rm
s}\times\ldots\times N_{N_{\rm D}}^{\rm s}}$ defines how the measurement of
the sparse signal is performed. To recover the coefficients tensor we can use
MOMP to solve:

$\min_{\bf C}\|{\bf y}-\sum_{{\bf i}\in\mathcal{I}}[{\bf\Phi}]_{:,{\bf
i}}\left[\prod_{k=1}^{N_{\rm
D}}{\bf\Psi}_{i_{k},j_{k}^{1}},\ldots,\prod_{k=1}^{N_{\rm
D}}{\bf\Psi}_{i_{k},j_{k}^{|\mathcal{J}|}}\right]{\bf C}_{\mathcal{J}}\|^{2}.$
(3)

Additionally, the support $\mathcal{C}\subset\mathcal{J}$ is defined as the
set of indices such that any element of ${\bf C}_{\mathcal{C}}$ is non zero.

###  III-B MOMP-based formulation

To formulate the mmWave channel estimation as a multidimensional orthogonal
matching pursuit problem, we need to rewrite (2) to match the structure of the
problem in (3). The next paragraphs derive the appropriate expressions for the
dictionaries and the sensing matrix so that MOMP can be used to solve (3)
given the observation matrix ${{\mathbf{Y}}}$.

The channel tensor ${\bf H}\in\mathbb{C}^{N_{\rm R}\times M_{\rm T}\times D}$
is defined as $[{\bf H}]_{:,:,d}={\bf H}_{d}$. To write the channel ${\bf H}$
as a sparse combination of products of dictionaries, it is natural to think of
independent dictionaries for the angle of arrival, the angle of departure, and
the delay, since these are the independent parameters of the paths that
compose the channel. Moreover, two different dictionaries are used to
represent the two independent dimensions of both the AoA and AoD. With this in
mind, we define five sparsifying dictionaries to represent the sparse mmWave
channel.

To obtain the compact expression of the dictionaries for the AoA, we consider
that for a horizontal uniform rectangular array of size $N_{\rm R}^{\rm
x}\times N_{\rm R}^{\rm y}$, the array response vector ${\bf a}_{\rm
R}(\bm{\theta}_{l})\in\mathbb{C}^{N_{\rm R}^{\rm x}N_{\rm R}^{\rm y}}$ can be
decomposed into two sub-components ${\bf a}_{\rm R}^{\rm x}(\theta^{\rm
x})\in\mathbb{C}^{N_{\rm R}^{\rm x}}$ and ${\bf a}_{\rm R}^{\rm y}(\theta^{\rm
y})\in\mathbb{C}^{N_{\rm R}^{\rm y}}$ such that ${\bf a}_{\rm
R}(\bm{\theta}_{l})={\bf a}_{\rm R}^{\rm x}(\theta^{\rm x})\otimes{\bf a}_{\rm
R}^{\rm y}(\theta^{\rm y})$. This means that we can rewrite the entries of
${\bf a}_{\rm R}(\bm{\theta}_{l})$ as $[{\bf a}_{\rm
R}(\bm{\theta}_{l})]_{n_{\rm R}^{\rm x}N_{\rm R}^{\rm y}+n_{\rm R}^{\rm
y}}=[{\bf a}_{\rm R}^{\rm x}(\theta_{l}^{\rm x})]_{n_{\rm R}^{\rm x}}[{\bf
a}_{\rm R}^{\rm y}(\theta_{l}^{\rm y})]_{n_{\rm R}^{\rm y}}$. Analogously, to
obtain the expression of the dictionaries for the AoD for a uniform
rectangular array of size $N_{\rm T}^{\rm x}\times N_{\rm T}^{\rm y}$ at the
transmitter, we can find the decomposition $[{\bf a}_{\rm
T}(\bm{\phi}_{l})]_{n_{\rm T}^{\rm x}N_{\rm T}^{\rm y}+n_{\rm T}^{\rm
y}}=[{\bf a}_{\rm T}^{\rm x}(\phi_{l}^{\rm x})]_{n_{\rm T}^{\rm x}}[{\bf
a}_{\rm T}^{\rm y}(\phi_{l}^{\rm y})]_{n_{\rm T}^{\rm y}}$. Finally, we will
define the dictionary for the delay domain from the evaluation of the pulse
shaping functions $[{\bf a}_{\rm D}(\tau)]_{d}=p((d-1)T_{\rm s}-\tau)$,
$d=1,\ldots,D$.

Now, to obtain the final expressions for the dictionaries, we can consider the
discrete domains for $\theta^{\rm x}$, $\theta^{\rm y}$, $\phi^{\rm x}$,
$\phi^{\rm y}$ and $\tau-\tau_{0}$ with resolutions $N^{\rm a}_{1}$, $N^{\rm
a}_{2}$, $N^{\rm a}_{3}$, $N^{\rm a}_{4}$ and $N^{\rm a}_{5}$, namely
$\\{\overline{\theta}_{1}^{\rm x},\ldots,\overline{\theta}_{N^{\rm
a}_{1}}^{\rm x}\\}$, $\\{\overline{\theta}_{1}^{\rm
y},\ldots,\overline{\theta}_{N^{\rm a}_{2}}^{\rm y}\\}$,
$\\{\overline{\phi}_{1}^{\rm x},\ldots,\overline{\phi}_{N^{\rm a}_{3}}^{\rm
x}\\}$, $\\{\overline{\phi}_{1}^{\rm y},\ldots,\overline{\phi}_{N^{\rm
a}_{4}}^{\rm y}\\}$ and $\\{\overline{\tau}_{1},\ldots,\overline{\tau}_{N^{\rm
a}_{5}}\\}$. This way, we define the dictionaries as

$\displaystyle{\bf\Psi}_{1}\\!=\\!\left[{\bf a}_{\rm R}^{\rm
x}(\overline{\theta}_{1}^{\rm x}),\ldots,{\bf a}_{\rm R}^{\rm
x}(\overline{\theta}_{N_{1}^{\rm e}}^{\rm
x})\right]\\!\\!,{\bf\Psi}_{2}\\!=\\!\left[{\bf a}_{\rm R}^{\rm
y}(\overline{\theta}_{1}^{\rm y}),\ldots,{\bf a}_{\rm R}^{\rm
y}(\overline{\theta}_{N_{2}^{\rm e}}^{\rm y})\right]\\!\\!,$
$\displaystyle{\bf\Psi}_{3}\\!\\!=\\!\\!\left[{\bf a}_{\rm T}^{\rm
x}(\overline{\theta}_{1}^{\rm x})^{*},\ldots,{\bf a}_{\rm T}^{\rm
x}(\overline{\theta}_{N_{3}^{\rm e}}^{\rm
x})^{*}\\!\right]\\!\\!,\\!{\bf\Psi}_{4}\\!\\!=\\!\\!\left[{\bf a}_{\rm
T}^{\rm y}(\overline{\theta}_{1}^{\rm y})^{*},\ldots,{\bf a}_{\rm T}^{\rm
y}(\overline{\theta}_{N_{4}^{\rm e}}^{\rm y})^{*}\\!\right]$
$\displaystyle{\bf\Psi}_{5}=\left[{\bf a}_{\rm
D}(\overline{\tau}_{1}),\ldots,{\bf a}_{\rm D}(\overline{\tau}_{N_{5}^{\rm
e}})\right].$ (4)

To use these dictionaries, we define the sets of multi-indicies
$\mathcal{I}=\\{{\bf i}=[i_{1},i_{2},i_{3},i_{4},i_{5}]\in\mathbb{N}\text{
such that }i_{1}\leq N_{\rm R}^{\rm x},i_{2}\leq N_{\rm R}^{\rm y},i_{3}\leq
N_{\rm T}^{\rm x},i_{4}\leq N_{\rm T}^{\rm y},i_{5}\leq D\\}$, and
$\mathcal{J}=\\{{\bf j}=[j_{1},j_{2},j_{3},j_{4},j_{5}]\in\mathbb{N}\text{
such that }j_{1}\leq N_{1}^{\rm a},j_{2}\leq N_{2}^{\rm a},j_{3}\leq
N_{3}^{\rm a},j_{4}\leq N_{4}^{\rm a},j_{5}\leq N_{5}^{\rm a}\\}$ and $N_{\rm
R}^{\rm x}\times N_{\rm R}^{\rm y}\times N_{\rm T}^{\rm x}\times N_{\rm
T}^{\rm y}\times D$ can be expressed as $\otimes_{k=1}^{5}N_{k}^{\rm s}$.
Finally, ignoring quantization effects caused by the finite resolution of the
dictionaries, we can define ${\bf
C}_{\mathcal{J}}\in\mathbb{C}^{\otimes_{d=1}^{5}N_{k}^{\rm e}\times 1}$ as the
sparse tensor

${\bf C}_{\mathcal{J}}=\left\\{\begin{array}[]{cl}\alpha_{l},&\text{if
}\begin{array}[]{ccc}\theta_{l}^{\rm x}=\overline{\theta}_{j_{1}}^{\rm x}\
{\rm or}\ \theta_{l}^{\rm y}=\overline{\theta}_{j_{2}}^{\rm y}\ {\rm or}\\\
\phi_{l}^{\rm x}=\overline{\phi}_{j_{3}}^{\rm x}\ {\rm or}\ \phi_{l}^{\rm
y}=\overline{\phi}_{j_{4}}^{\rm y}\ {\rm or}\
\tau_{l}=\overline{\tau}_{j_{5}}\end{array}\\\
0,&\text{otherwise.}\end{array}\right.$ (5)

Under these definitions, the channel entries can now be expressed as

$[{\bf H}]_{i_{1}N_{\rm R}^{\rm y}+i_{2},i_{3}N_{\rm T}^{\rm
y}+i_{4},i_{5}}=\left[\prod_{k=1}^{N_{\rm
D}}{\bf\Psi}_{i_{k},j_{k}^{1}},\ldots,\prod_{k=1}^{N_{\rm
D}}{\bf\Psi}_{i_{k},j_{k}^{|\mathcal{J}|}}\right]{\bf C}_{\mathcal{J}}.$ (6)

Regarding the construction of the observation, we assume that each training
frame is independent. Next, the received signal is whitened to compensate for
the correlation effect introduced by the analog combiner. To this aim, we find
a Cholesky decomposition of the noise correlation matrix ${\bf W}_{m}^{\rm
H}{\bf W}_{m}$ i.e. ${\bf L}_{m}$ such that ${\bf L}_{m}{\bf L}_{m}^{\rm
H}={\bf W}_{m}^{\rm H}{\bf W}_{m}$ and multiply (2) by ${\bf L}_{m}^{-1}$ to
obtain

$[{\bf L}_{m}^{-1}{\bf Y}_{m}]_{:,q}=\sqrt{P}\sum_{d=1}^{D}{\bf
L}_{m}^{-1}{\bf W}_{m}^{\rm H}[{\bf H}]_{:,:,d}{\bf F}_{m}[{\bf
S}_{m}]_{:,q+D-d}\\\ +{\bf L}_{m}^{-1}{\bf W}_{m}^{\rm H}[{\bf N}_{m}]_{:,q}.$
(7)

After the whitening stage, the expressions for the observation ${\bf y}$, and
noise ${\bf n}\in\mathbb{C}^{MM_{\rm R}Q}$ are

$\displaystyle[{\bf y}]_{mM_{\rm R}Q+m_{\rm R}Q+q}=$ $\displaystyle[{\bf
L}_{m}^{-1}{\bf Y}_{m}]_{m_{\rm R},q}$ (8) $\displaystyle[{\bf n}]_{mM_{\rm
R}Q+m_{\rm R}Q+q}=$ $\displaystyle[{\bf L}_{m}^{-1}{\bf W}_{m}^{\rm H}{\bf
N}_{m}]_{m_{\rm R},q}.$ (9)

To complete the MOMP formulation, we still need to define the sensing matrix
${\bf\Phi}\in\mathbb{C}^{MM_{\rm R}Q\times N_{\rm R}^{\rm x}\times N_{\rm
R}^{\rm y}\times N_{\rm T}^{\rm x}\times N_{\rm T}^{\rm y}\times D}$. Note
that we can also write the observations in (8) as $[{\bf y}]_{mM_{\rm
R}Q+m_{\rm R}Q+q}=[{\bf\Phi}]_{mM_{\rm R}Q+m_{\rm R}Q+q,{\bf i}}[{\bf
H}]_{i_{1}N_{\rm R}^{\rm y}+i_{2},i_{3}N_{\rm T}^{\rm y}+i_{4},i_{5}}$.
Mapping the terms in this observation to (7)

$[{\bf\Phi}]_{mM_{\rm R}Q+m_{\rm R}Q+q,{\bf i}}=\\\ \sqrt{P}[{\bf
L}_{m}^{-1}{\bf W}_{m}^{\rm H}]_{m_{\rm R},i_{1}N_{\rm R}^{\rm y}+i_{2}}[{\bf
F}_{m}{\bf S}_{m}]_{i_{3}N_{\rm T}^{\rm y}+i_{4},q+D-i_{5}}.$ (10)

With these definitions, the multi-dimensional matching pursuit problem is
completed, since we can re-write equation (2) like

${\bf y}=\sum_{{\bf i}\in\mathcal{I}}[{\bf\Phi}]_{:,{\bf
i}}\left[\prod_{k=1}^{N_{\rm
D}}{\bf\Psi}_{i_{k},j_{k}^{1}},\ldots,\prod_{k=1}^{N_{\rm
D}}{\bf\Psi}_{i_{k},j_{k}^{|\mathcal{J}|}}\right]{\bf C}_{\mathcal{J}}+{\bf
n}.$ (11)

Since ${\bf n}$ is white noise, the maximum likelihood estimator is the
solution to (3), and we can apply the MOMP algorithm to approximately and
sparsely solve it.

The parameters for the estimated path corresponding to the index ${\bf
j}\in\mathcal{C}$ can be directly extracted from the support and the sparse
reconstructed matrix as $\overline{\theta}_{j_{1}}^{\rm x}$,
$\overline{\theta}_{j_{2}}^{\rm y}$, $\overline{\phi}_{j_{3}}^{\rm x}$,
$\overline{\phi}_{j_{4}}^{\rm y}$ and $\overline{\tau}_{j_{5}}$ for the
angular information and relative time of arrival, and as ${\bf
C}_{\mathcal{J}}$ for the complex gain. To fully retrieve the angular
information, the z components of $\bm{\theta}$ and $\bm{\phi}$ can be computed
using the fact that these are unitary vectors, and that the z component is
positive (otherwise the path would be coming from the antenna array substrate
which blocks the signal), i.e. $\overline{\theta}^{\rm
z}=\sqrt{(\overline{\theta}^{\rm x})^{2}+(\overline{\theta}^{\rm y})^{2}}$ and
$\overline{\phi}^{\rm z}=\sqrt{(\overline{\phi}^{\rm
x})^{2}+(\overline{\phi}^{\rm y})^{2}}$.

##  IV Localization

Once the channel parameters have been estimated, the geometric relationships
between those parameters and the scatterers in the environment can be used to
obtain an estimate of the transmitter position. These geometric relationships
between the channel parameters and the position to be estimated depend on the
type of path. Therefore, a method to classify the different estimated paths
has to be proposed. To this aim, we will make use of two properties: a)
specular reflections on the horizontal plane (floor/ceiling) do not change the
$x$-$y$ components of a point, and b) vertical reflections (walls) do not
change the $z$ component. With this in mind, we consider 4 possible types of
paths: line of sight (LoS), wall reflection, floor/ceiling reflection, or any
other path that will be labeled as spurious and will not be exploited for
localization.

The input to our proposed path classification algorithm are the spherical
coordinates of the different angles, computed as $\theta^{\rm
az}=\arg(\theta^{\rm x}+j\theta^{\rm y})$, $\theta^{\rm
el}=\arcsin(\theta^{\rm z})$, $\phi^{\rm az}=\arg(\phi^{\rm x}+j\phi^{\rm y})$
and $\phi^{\rm el}=\arcsin(\phi^{\rm z})$. Wall reflections and LoS paths
satisfy $\theta_{l}^{\rm el}+\phi_{l}^{\rm el}=0$ while floor/ceiling
reflections satisfy $\theta_{l}^{\rm el}=\phi_{l}^{\rm el}$. Additionally,
floor/ceiling reflections and LoS paths arrival and departure azimuth angles
are opposite, i.e. they are separated by $180^{\circ}$. By defining the
threshold values $r_{\rm az},r_{\rm el}$ we can obtain the conditions

$\displaystyle|\sin(\theta^{\rm el}-\phi_{l}^{\rm el})|<r_{\rm el},$ (12)
$\displaystyle|\sin(\theta^{\rm el}+\phi_{l}^{\rm el})|<r_{\rm el},$ (13)
$\displaystyle\cos(\theta^{\rm az}-\phi_{l}^{\rm az})<r_{\rm az}-1.$ (14)

A LoS path satiifies conditions (12) and (14), floor/ceiling reflections
satisfy (13) and (14), while wall reflections only satisfy (12). Any other
path will be classified as spurious.

Let us define now the algorithm that provides the position estimate from the
path parameters and path classification. We consider the access point
(receiver) to be in the coordinate origin, and ${\bf u}\in\mathbb{R}^{3}$ to
be the device (transmitter) location. To ease the formulation of the position
estimation algorithm we introduce the projection matrices ${\bm{\chi}}_{\rm
LOS}={\bf I}_{3}$, ${\bm{\chi}}_{\rm h}=\left[\begin{array}[]{cc}1&0\\\ 0&1\\\
0&0\end{array}\right]$, and ${\bm{\chi}}_{\rm v}=\left[\begin{array}[]{c}0\\\
0\\\ 1\end{array}\right]$. With these definitions, the geometric relationships
created by a LoS path, a floor/ceiling reflection or a wall reflection can be
expressed as

$\displaystyle{\bm{\chi}}_{\rm LOS}{\bf u}=\\!{\bm{\chi}}_{\rm
LOS}{\bm{\theta}}_{l}\tau_{l},\ {\bm{\chi}}_{\rm h}{\bf u}=\\!{\bm{\chi}}_{\rm
h}{\bm{\theta}}_{l}\tau_{l},\ {\bm{\chi}}_{\rm v}{\bf u}=\\!{\bm{\chi}}_{\rm
v}{\bm{\theta}}_{l}\tau_{l},$ (15)

respectively. In other words, if we define ${\bm{\chi}}_{l}$, as the
projection matrix corresponding to the $l$-th path type, we can rewrite the
geometric relationships for that path as

${\bm{\chi}}_{l}({\bf u}-{\bm{\theta}}_{l}\tau_{l})=0.$ (16)

Note, however, that the channel estimation algorithm only provides the
estimation of the relative delays $\Delta\tau_{l}=\tau_{l}-\tau_{0}$, with
$\tau_{0}$ unknown, instead of $\tau_{l}$. Therefore, the equations for the
geometric relationships need to be established in terms of the relative
delays, i.e.

${\bm{\chi}}_{l}({\bf u}-{\bm{\theta}}_{l}(\Delta\tau_{l}+\tau_{0}))=0.$ (17)

Now, we define the matrix ${\bf U}_{{\bm{\theta}}_{l}}=[{\bf
I},-{\bm{\theta}}_{l}]$ and the variable ${\bf z}=[{\bf u}^{\rm
T},\tau_{0}]^{\rm T}$, which contains all the unknown variables, to rewrite
(17) as

${\bm{\chi}}_{\rm p}({\bf U}_{{\bm{\theta}}_{l}}{\bf
z}-{\bm{\theta}}_{l}\Delta\tau_{l})=0.$ (18)

This equation can be alternatively solved by

$\min_{{\bf z}}{\bf z}^{\rm T}{\bf U}_{{\bm{\theta}}_{l}}^{\rm
T}{\bm{\chi}}_{l}^{\rm T}{\bm{\chi}}_{l}{\bf U}_{{\bm{\theta}}_{l}}{\bf
z}-2\tau_{l}\phi_{l}^{\rm T}{\bm{\chi}}_{l}^{\rm T}{\bm{\chi}}_{l}{\bf
U}_{{\bm{\theta}}_{l}}{\bf z}+\tau_{l}^{2}\phi_{l}^{\rm T}{\bm{\chi}}_{l}^{\rm
T}{\bm{\chi}}_{l}\phi_{l}.$ (19)

We define now a weighted version of the previous equation to account for the
different impact of the errors associated to different types of paths. Thus,
we consider a weight value for each path $w_{l}$ that we define as its
estimated gain $w_{l}=|\alpha_{l}|^{2}$. With this in mind we can weight
average (19) for all the different paths as

$\min_{{\bf z}}{\bf z}^{\rm T}{\bf A}{\bf z}-2{\bf b}{\bf z}+c,$ (20)

for ${\bf A}=\sum_{l=1}^{L}w_{l}{\bf U}_{{\bm{\theta}}_{l}}^{\rm
T}{\bm{\chi}}_{l}^{\rm T}{\bm{\chi}}_{l}{\bf U}_{{\bm{\theta}}_{l}}$, ${\bf
b}=\sum_{l=1}^{L}w_{l}\tau_{l}\phi_{l}^{\rm T}{\bm{\chi}}_{l}^{\rm
T}{\bm{\chi}}_{l}{\bf U}_{{\bm{\theta}}_{l}}$ and
$c=\sum_{l=1}^{L}w_{l}\tau_{l}^{2}\phi_{l}^{\rm T}{\bm{\chi}}_{l}^{\rm
T}{\bm{\chi}}_{l}\phi_{l}$. The widely known solution to (20) for a symmetric
${\bf A}$ matrix is $\hat{\bf z}={\bf A}^{-1}{\bf b}$. The location estimation
$\hat{\bf u}$ can then be extracted from the first 3 entries of $\hat{\bf z}$,
while the last entry provides the unkown offset between transmitter and
receiver.

##  V Simulation

To generate the evaluation data set with channels and user positions we
simulate by ray tracing (using Wireless InSite sofware) a home office with two
rooms. There is an access point in each one of the rooms. We generate $218$
user locations following a path connecting both rooms. The user antenna is
place horizontally, while the access points are mounted on the walls
vertically, facing the interior of their respective rooms. The ray tracing
software generates the optical geometrical channel between the user positions
and the access points. We associate each user location to the single access
point with the highest received signal strength. For each user location, we
input the parameters corresponding to the channel with its associated access
point in a measurement model to generate the received signal. The transmit
power is set to $P=20{\rm dBm}$, and the noise level is the thermal noise
corresponding to a room at $15^{\circ}{\rm C}$ when using a $1{\rm GHz}$
bandwidth, i.e. $\sigma^{2}=-84{\rm dBm}$. The pulse shaping filter is defined
as a sinc function. We set the delay tap length to $D=64$, and the training
signal to be the first $6$ rows of a $64$ element Hadamard matrix padded with
$64$ zeros to the left and $32$ zeros to the right, i.e. ${\bf
S}=[0,\ldots,0,{\bf Had}_{64},0,\ldots,0]$. The training precoders/combiners
are taken as columns extracted from the Kronecker product of the DFT matrices
with sizes $N_{\rm T}^{\rm x}$ and $N_{\rm T}^{\rm y}$/$N_{\rm R}^{\rm x}$ and
$N_{\rm R}^{\rm y}$.

To simplify the sparse reconstruction problem, we define the ratio between the
number of atoms and atoms size $K_{\rm res}=N_{k}^{\rm a}/N_{k}^{\rm s}$,
setting it to the same constant value independent of $k$. By how the
dictionaries are defined we have $N_{1}^{\rm s}=N_{\rm R}^{\rm x}$,
$N_{2}^{\rm s}=N_{\rm R}^{\rm y}$, $N_{3}^{\rm s}=N_{\rm T}^{\rm x}$,
$N_{4}^{\rm s}=N_{\rm R}^{\rm y}$, and $N_{5}^{\rm s}=D$.

The first metric used to evaluate the accuracy of the channel estimation
algorithm is the error in the estimation of the AoA of the main path. We can
see its evolution over $K_{\rm res}$ for both algorithms, OMP and MOMP, in
Table I, when considering a $3\times 3$ antenna array at the transmitter and a
$6\times 6$ antenna at the receiver. The number of RF chains is set to 3 for
the device and to 6 for the access point.

Method | OMP | MOMP
---|---|---
$K_{\rm res}$ | $1$ | $1.3$ | $1.6$ | $16$ | $128$ | $1024$
$\angle(\hat{\bm{\theta}}_{1},{\bm{\theta}}_{1})[^{\circ}]$ | $12$ | $10$ | $8$ | $1.1$ | $0.9$ | $0.9$
$\angle(\hat{\bm{\phi}}_{1},{\bm{\phi}}_{1})[^{\circ}]$ | $44$ | $30$ | $25$ | $2.1$ | $1.8$ | $1.7$
Run-time $[{\rm s}]$ | $1.3$ | $2.2$ | $8.2$ | $3.2$ | $3.3$ | $3.5$
TABLE I: Median main path angular error and run time as a function of $K_{\rm
res}$.

It is straight forwards to see how MOMP outperforms OMP estimations by an
order of magnitude because of its ability to increase the resolution of the
dictionaries without incurring in exorbitant computational times or impossible
memory requirements.

We analyze next the localization error $\|{\bf u}-\hat{\bf u}\|$ as a function
of the number of antennas and complexity. The value used for the weight is
each path’s estimated power $w_{l}=|\bar{\alpha}_{l}|^{2}$ and the threshold
values for the path classification are $r_{\rm az}=r_{\rm el}=0.12$. We fix
the antenna size for the device to $3\times 3$, while three different sizes
are considered for the access points: $6\times 6$, $8\times 8$, and $10\times
10$. The number of receive RF chains also varies accordingly as 6, 8 and 10.
The localization results can be visualized in Fig. 1. Even with these small
antennas sizes, the access point is able to localize most of the devices
keeping the error under reasonable values when an unknown clock offset is
considered in the simulation. Note that even channels with very bad SNR
conditions, where the link cannot be established have been considered in these
averaged results. These results clearly show the benefits of applying MOMP to
the problem of joint channel estimation and localization at mmWave. Higher
estimation accuracies are expected with larger antenna arrays at the
transmitter, which is constrained to be only $3\times 3$ in this simulation
scenario so OMP can run.

Figure 1: Localization error as a function of the receive array size, number
of RF chains and sparse recovery strategy. $K_{\rm res}$ has been set to 128
for MOMP and to 1.6 for OMP.

##  VI Conclusions

We developed a joint compressive channel estimation and localization strategy
for a realistic mmWave MIMO communication systems. The proposed approach is
based on a multidimensional matching pursuit sparse recovery algorithm, which
enables operation with high resolution dictionaries, for both the time and
angular domains, at a reduced complexity. Using mmWave channels generated by
ray tracing we obtained the channel estimation error for the angular
parameters and the localization errors for a random deployment of users, as a
function of the complexity and the system parameters. We showed how the MOMP-
based approach is able to provide reasonable estimation accuracies in
scenarios where conventional OMP is not feasible.

## References

  * [1] A. Shahmansoori, G. E. Garcia, G. Destino, G. Seco-Granados, and H. Wymeersch, “Position and orientation estimation through millimeter-wave MIMO in 5G systems,” _IEEE Trans. on Wireless Commun/_ , vol. 17, no. 3, pp. 1822–1835, 2018.
  * [2] H. Wymeersch, N. Garcia, H. Kim, G. Seco-Granados, S. Kim, F. Went, and M. Fröhle, “5G mmWave downlink vehicular positioning,” in _IEEE Global Commun. Conf. (GLOBECOM)_. IEEE, 2018, pp. 206–212.
  * [3] J. Talvitie, M. Koivisto, T. Levanen, M. Valkama, G. Destino, and H. Wymeersch, “High-accuracy joint position and orientation estimation in sparse 5G mmWave channel,” in _IEEE Intl. Conf. on Commun. (ICC)_ , 2019, pp. 1–7.
  * [4] F. Jiang, Y. Ge, M. Zhu, and H. Wymeersch, “High-dimensional channel estimation for simultaneous localization and communications,” in _2021 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2021, pp. 1–6.
  * [5] K. Venugopal, A. Alkhateeb, N. Prelcic-González, and R. W. Heath, “Channel estimation for hybrid architecture-based wideband millimeter wave systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 35, no. 9, pp. 1996–2009, 2017.
  * [6] J. Rodríguez-Fernández, N. González-Prelcic, K. Venugopal, and R. W. Heath, “Frequency-domain compressive channel estimation for frequency-selective hybrid millimeter wave MIMO systems,” _IEEE Trans. on Wireless Commun._ , vol. 17, no. 5, pp. 2946–2960, 2018.
  * [7] X. Wu, G. Yang, F. Hou, and S. Ma, “Low-complexity downlink channel estimation for millimeter-wave FDD massive MIMO systems,” _IEEE Wireless Communications Letters_ , 2019.
  * [8] F. Zhu, A. Liu, and V. Lau, “Channel estimation and localization for mmWave systems: A sparse bayesian learning approach,” in _IEEE Intl. Conf. on Commun. (ICC)_ , 2019, pp. 1–6.
  * [9] R. W. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 10, no. 3, pp. 436–453, April 2016.
  * [10] J. Palacios, N. González-Prelcic, and C. Rusu, “Multidimensional orthogonal matching pursuit: theory and application to joint channel estimation and localization at mmWave,” arXiv preprint, 2022.
  * [11] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” _IEEE Transactions on Information Theory_ , vol. 53, no. 12, pp. 4655–4666, 2007.