Power Allocation and Parameter Estimation for Multipath-based 5G Positioning
Anastasios Kakkavas, Henk Wymeersch, Gonzalo Seco-Granados, Mario H. Castañeda García, Richard A. Stirling-Gallacher, Josef A. Nossek
11 Power Allocation and Parameter Estimation forMultipath-based 5G Positioning
Anastasios Kakkavas,
Student Member, IEEE,
Henk Wymeersch,
Senior Member, IEEE,
Gonzalo Seco-Granados,
Senior Member, IEEE,
Mario H. Castañeda García,
Member, IEEE,
Richard A. Stirling-Gallacher,
Member, IEEE, and Josef A. Nossek,
Life Fellow, IEEE
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible.
Abstract
We consider a single-anchor multiple-input multiple-output orthogonal frequency-division multi-plexing system with imperfectly synchronized transmitter (Tx) and receiver (Rx) clocks, where theRx estimates its position based on the received reference signals. The Tx, having (imperfect) priorknowledge about the Rx location and the surrounding geometry, transmits the reference signals basedon a set of fixed beams. In this work, we develop strategies for the power allocation among the beams
This work was supported in part by the EU-H2020 project Fifth Generation Communication Automotive Research andInnovation (5GCAR), and in part by the ICREA Academia program and the Spanish Ministry of Science, Innovation andUniversities project TEC2017-89925-R.A. Kakkavas is with the Munich Research Center, Huawei Technologies Duesseldorf GmbH, 80992 Munich, Germany, andalso with the Department of Electrical and Computer Engineering, Technical University of Munich, 80333 Munich, Germany(e-mail: [email protected]).H. Wymeersch is with the Department of Electrical Engineering, Chalmers University of Technology, 412 58 Gothenburg,Sweden (email: [email protected]).G. Seco-Granados is with the Department of Telecommunications and Systems Engineering, Universitat Autonoma deBarcelona, Spain (UAB) (e-mail: [email protected]).M. H. Castañeda García and R. A. Stirling-Gallacher are with the Munich Research Center, Huawei Technologies DuesseldorfGmbH, 80992 Munich, Germany (e-mail: [email protected]; [email protected]).J. A. Nossek is with the Department of Electrical and Computer Engineering, Technical University of Munich, 80333 Munich,Germany (e-mail: [email protected]). a r X i v : . [ c s . I T ] D ec aiming to minimize the expected Cramér-Rao lower bound for Rx positioning. Additional constraintson the design are included to ensure that the line-of-sight (LOS) path is detected with high probability.Furthermore, the effect of clock asynchronism on the resulting allocation strategies is also studied.We also propose a gridless compressed sensing-based position estimation algorithm, which exploits theinformation on the clock offset provided by non-line-of-sight paths, and show that it is asymptoticallyefficient. Index Terms positioning, localization, 5G, reference signal, power allocation, parameter estimation
I. I
NTRODUCTION
With the advent of fifth generation (5G) mobile networks, positioning has attracted lotsof research interest. The large chunks of bandwidth available at millimeter-wave (mm-Wave)frequencies, as well as the potentially large number of antennas placed at both sides of thecommunication link are the main driving forces, not only for very high data rates and massiveconnectivity [1], [2], but also for a drastic improvement of the positioning accuracy of cellularnetworks [3]. Recently, within the Third Generation Partnership Project (3GPP), besides po-sitioning techniques already existing in previous generations of cellular networks [4], such asobserved time difference of arrival (OTDOA), uplink TDOA (UTDOA), new techniques havebeen standardized, including downlink (DL)-angle of departure (AOD), uplink (UL)-angle ofarrival (AOA) and multi-cell round-trip time (RTT) [5]. In addition, proposals for reportingdelay and angular multipath measurements to enable single-anchor positioning have been con-sidered [6]. With their enhanced positioning capabilities, 5G systems aim to accommodate usecases like assisted/autonomous driving [7], augmented reality and industrial internet of things(IoT) [6].Single-anchor localization, that leverages the high temporal and angular resolution of mm-Wave multiple-input multiple-output (MIMO) systems, has received increasing attention in recentyears, as it has the potential to ease the requirements of multi-anchor hearability and interferencemanagement. The fundamental limits of single-anchor positioning have been investigated in [8]for line-of-sight (LOS) and [9]–[12] for multipath channels with single-bounce-non-LOS (NLOS)components.The single-anchor localization algorithms in the literature can be classified into two categories:one-shot schemes without tracking [13]–[20], and tracking approaches [21]–[27]. While the latter mainly focus on positon estimation and tracking given the channel parameter measurements, theformer also deal with the estimation of the channel parameters, as done in the present work.A three-stage algorithm for position estimation with a multiple-input multiple-output (MIMO)-orthogonal frequency-division multiplexing (OFDM) system was proposed in [13], where inthe first stage a compressed sensing-based algorithm is used to obtain coarse estimates of themultipath parameters (number of paths, times of arrival (TOAs), AODs, AOAs and gains), withthe coarse estimates refined in the second stage. In the third stage, the refined estimates aremapped to the receiver (Rx) position and orientation and the scatterer/reflector positions usingthe extended invariance principle (EXIP). A similar approach is followed in [14], with themain difference lying in the mapping from channel parameters to position parameters, wherean iterative Gibbs sampling method is employed. In [15] range-free angle-based approaches aredeveloped assuming prior map information. An algorithm for localization and synchronizationof cooperating full-duplex agents using a single-anchor is developed in [16]. The authors of [17]propose a protocol and an accompanying algorithm that enables a single-anchor to (quasi-)simultaneously receive messages from multiple agents in order to localize them using TOAand AOA measurements. The proposed approach is verified on an experimental setup. A DLpositioning algorithm for a single-antenna Rx, based on TOA and AOD measurements is proposedin [19]. The work is extended in [20], where a two-step process is used, with the coarse parameterestimates obtained in the first step used for adaptation of the transmitter (Tx) beamforming matrixin the second step. In [28] an iterative position estimation and Tx beamforming refinementalgorithm is developed.Similar to [20], [28], many works have considered the use of prior knowledge of the Rxposition at the Tx to design beamformers that improve the Rx’s localization accuracy. In [29]Cramér-Rao lower bound (CRLB)-optimal precoders for tracking the AOD and AOA of a pathwere designed, taking the uncertainty about their value into account. In [30], assuming a LOSchannel and a multicarrier system, beamformers minimizing the TOA and AOA error boundswere proposed, based on the current estimate of the Rx position. Using a similar setup, butadditionally considering multiple users, the authors of [31] designed beamformers maximizing aweighted sum of Fisher information on delay, AOD and AOA. Although in a different context,the algorithms and the conclusions of [32], [33] are relevant to our Tx beamforming problem.In [32], [33], robust beamformers under angular uncertainty were designed and it is concludedthat the Rx steering vector and its derivative contain all the localization information. Again in a different but still relevant setup, the authors of [34] and [35] compute the optimal power allocationamong multiple anchors for ranging-based localization by solving a semidefinite program (SDP).In [36] it was shown that, when the uncertainty about the Rx position is not considered, it isoptimal to transmit only on the directions corresponding to the Tx array steering vector and itsderivative. The power allocation among these two directions minimizing the squared positionerror bound (SPEB) qas analytically calculated in [36]. When the Rx location uncertainty is takeninto account, the optimal power allocation among the beams of a given Tx beam codebook wascomputed to minimize the average or maximum SPEB.In this paper, we extend our work in [36]. We consider a single-anchor setup and a sparsemultipath channel, which comprises the LOS path and a number of single-bounce NLOS paths,as multi-bounce paths are considered too weak for reception at mm-Wave frequencies [37]–[40].The Tx has only a coarse prior knowledge of the underlying geometry and in addition, the Tx-Rx clocks are imperfectly synchronized. We optimize the power allocation on a beam codebookfor the multipath channel and examine the effect of imperfect synchronization on the resultingpower allocation. Also, we develop a novel position estimation algorithm, which is evaluated forthe proposed power allocation strategies. The main contributions of the work can be summarizedas follows: • We propose power allocation strategies on a fixed Tx beam codebook with the aim ofminimizing the expected positioning error of the Rx. The optimal solution and a suboptimalone with lower computational complexity are presented and evaluated. • We develop a two-stage position estimation algorithm. The first stage consists of a gridlesschannel parameter estimation algorithm, based on [41]. The second stage maps the channelparameter estimates to position parameters. The information about the clock offset offeredby NLOS paths in combination with the LOS path is exploited so as to discard false alarms.The rest of the paper is organized as follows. In Sec. II we present the system modeland the assumptions of the work. The theoretical bound on positioning accuracy is brieflydiscussed in Sec. III and the proposed power allocation methods are presented in Sec. IV. Theposition estimation algorithm is introduced in Sec. V and numerical evaluations of the proposedapproaches are provided in Sec. VI. Sec. VII concludes the work.
Notation:
We use bold lowercase for vectors, bold uppercase for matrices, non-bold for scalarsand calligraphic letters for sets. Depending on its argument, | · | denotes the absolute value of ascalar, the determinant of a matrix or the cardinality of a set. The transpose, conjugate transpose p T p T , j ψ T , j d T , j p R α R p R , i ψ R , i d R , i θ T , p s , θ R , p s , Fig. 1. Geometric model, example with a uniform linear array (ULA) at the Tx and a uniform circular array (UCA) at the Rx. and 𝑝 -norm of a vector/matrix are denoted by (·) T , (·) H and (cid:107) · (cid:107) 𝑝 and the Frobenius norm ofa matrix is denoted by (cid:107) · (cid:107) F . (cid:60){·} and (cid:61){·} denote the real and imaginary part of a complexnumber and arg (·) denotes its phase. The 𝑖 -th element of a vector and the ( 𝑖, 𝑗 ) -th element of amatrix are denoted by [·] 𝑖 and [·] 𝑖, 𝑗 , respectively. I 𝑛 , and denote the identity matrix of size 𝑛 , and the all-ones all-zeros matrix of the appropriate size. diag ( x ) denotes the diagonal matrixwith the elements of x on its diagonal. The expectation operator is denoted by E [·] and thesets of real and complex numbers are denoted by R and C . A multivariate (circularly symmetriccomplex) Gaussian distribution with mean µ and covariance matrix C is denoted by N ( µ , C ) ( N C ( µ , C ) ). The Hessian of a function 𝑓 ( x ) is denoted as 𝐷 x 𝑓 ( x ) .II. S YSTEM M ODEL AND A SSUMPTIONS
A. Geometric Model
The Tx consists of an array with 𝑁 T antennas and reference point located at the origin. TheRx consists of an array with 𝑁 R antennas, a reference point located at p R = (cid:2) 𝑝 R,x , 𝑝
R,y (cid:3) T ∈ R and orientation 𝛼 R . The position of the 𝑗 -th element of the Tx array is given by p T , 𝑗 = 𝑑 T , 𝑗 u (cid:0) 𝜓 T , 𝑗 (cid:1) ∈ R , 𝑗 = , . . . , 𝑁 T − , (1)where u (cid:0) 𝜓 (cid:1) = (cid:2) cos (cid:0) 𝜓 (cid:1) , sin (cid:0) 𝜓 (cid:1) (cid:3) T and 𝑑 T , 𝑗 and 𝜓 T , 𝑗 are its distance and angle from the Tx array’sreference point, as shown in Fig. 1. Accordingly, the position of the 𝑖 -th element of the Rx arrayis given by p R ,𝑖 = 𝑑 R ,𝑖 u ( 𝜓 R ,𝑖 + 𝛼 R ) ∈ R , 𝑖 = , . . . , 𝑁 R − . (2)We assume that for all antenna pairs there are 𝐿 discrete propagation paths. The first of these 𝐿 paths ( 𝑙 = ) is the LOS path and the rest ( 𝑙 = , . . . , 𝐿 − ) are single-bounce NLOS paths. The point of incidence of the 𝑙 -th single-bounce path, which corresponds either to scatteringor reflection, is p s ,𝑙 = (cid:2) 𝑝 s ,𝑙, x , 𝑝 s ,𝑙, y (cid:3) T , 𝑙 = , . . . , 𝐿 − . The array apertures are assumed to besmall compared to the distance between Tx and Rx, as well as the distance between each of thescatterers/reflectors and the Tx or Rx. Therefore, the delay of the 𝑙 -th path from Tx element 𝑗 to Rx element 𝑖 can be approximated by [12] 𝜏 𝑙,𝑖, 𝑗 ≈ 𝜏 𝑙 − 𝜏 T , 𝑗 ( 𝜃 T ,𝑙 ) − 𝜏 R ,𝑖 ( 𝜃 R ,𝑙 ) , 𝑙 = , . . . , 𝐿 − , (3)where 𝜏 𝑙 = (cid:107) p R (cid:107) / 𝑐 + 𝜖 clk , 𝑙 = (cid:0) (cid:107) p s ,𝑙 (cid:107) + (cid:107) p R − p s ,𝑙 (cid:107) (cid:1) / 𝑐 + 𝜖 clk , 𝑙 ≠ , (4) 𝜏 T , 𝑗 ( 𝜃 T ,𝑙 ) = 𝑑 T , 𝑗 u T ( 𝜓 T , 𝑗 ) u ( 𝜃 T ,𝑙 )/ 𝑐, (5) 𝜏 R ,𝑖 ( 𝜃 R ,𝑙 ) = 𝑑 R ,𝑖 u T ( 𝜓 R ,𝑖 ) u ( 𝜃 R ,𝑙 )/ 𝑐, (6)with 𝜖 clk being the clock offset between Tx and Rx and 𝑐 the speed of light. The angles aredefined as 𝜃 T ,𝑙 = atan2 (cid:16) 𝑝 R , y , 𝑝 R , x (cid:17) , 𝑙 = (cid:16) 𝑝 s ,𝑙, y , 𝑝 s ,𝑙, x (cid:17) , 𝑙 ≠ (7) 𝜃 R ,𝑙 = 𝜃 T ,𝑙 + 𝜋 − 𝛼 R , 𝑙 = , atan2 (cid:16) 𝑝 s ,𝑙, y − 𝑝 R , y , 𝑝 s ,𝑙, x − 𝑝 R , x (cid:17) − 𝛼 R , 𝑙 ≠ , (8)with atan2 (cid:0) 𝑦, 𝑥 (cid:1) being the four-quadrant inverse tangent function. B. Signal Model
An OFDM waveform with subcarrier spacing Δ 𝑓 , 𝑁 subcarriers and cyclic prefix (CP)duration 𝑇 CP is considered. The reference signal is transmitted on 𝑁 P subcarriers, whose indicesare described by P = { 𝑝 , . . . , 𝑝 𝑁 P } and 𝑁 B OFDM symbols are transmitted. We assume anarrowband signal model, i.e. 𝐵 / 𝑓 c (cid:28) 𝜆 c / 𝐷 max , where 𝐵 ≈ Δ 𝑓 ( max (P) − min (P)) is the signalbandwidth, 𝑓 𝑐 is the carrier frequency, 𝜆 c is the carrier wavelength and 𝐷 max is the largest of theTx and Rx array apertures. The reference signal resource grid R comprises all resource elementsat the time-frequency points ( 𝑝, 𝑏 ) , 𝑝 ∈ P , 𝑏 = , . . . , 𝑁 B − . The transmitter uses a beamcodebook { f 𝑘 } 𝑀 T 𝑘 = , where 𝑀 T is the number of beams in the codebook and (cid:107) f 𝑘 (cid:107) = , ∀ 𝑘 . The 𝑘 -th beam is used on a subset R 𝑘 of resource elements (REs) ( 𝑝, 𝑏 ) , with R 𝑘 ∩ R 𝑘 (cid:48) = ∅ for 𝑘 ≠ 𝑘 (cid:48) . The transmitted signal vector at the 𝑝 -th subcarrier, 𝑝 ∈ P , of the 𝑏 -th OFDM symbol, 𝑏 = , . . . , 𝑁 B − , then is x [ 𝑝, 𝑏 ] = 𝜆 𝑘 [ 𝑝, 𝑏 ] f 𝑘 , ( 𝑝, 𝑏 ) ∈ R 𝑘 , (9)where 𝜆 𝑘 [ 𝑝, 𝑏 ] = √︁ 𝑃 tot 𝑞 𝑘 𝛾 𝑘 [ 𝑝, 𝑏 ] e j 𝛽 𝑘 [ 𝑝,𝑏 ] (10)is the symbol assigned to f 𝑘 at the 𝑝 -th subcarrier, 𝑃 tot is the total Tx power (disregarding thepower spent for the CP), 𝑞 𝑘 is the fraction of 𝑃 tot allocated to f 𝑘 , with (cid:205) 𝑀 T 𝑘 = 𝑞 𝑘 = , 𝛾 𝑘 [ 𝑝, 𝑏 ] is the fraction of 𝑞 𝑘 allocated to the RE ( 𝑝, 𝑏 ) , with (cid:205) ( 𝑝,𝑏 )∈R 𝑘 𝛾 𝑘 [ 𝑝, 𝑏 ] = , and 𝛽 𝑘 [ 𝑝, 𝑏 ] is thephase of 𝜆 𝑘 [ 𝑝, 𝑏 ] . The received signal is y [ 𝑝, 𝑏 ] = m [ 𝑝, 𝑏 ] + η [ 𝑝, 𝑏 ] , (11)where m [ 𝑝, 𝑏 ] = 𝐿 − ∑︁ 𝑙 = ℎ 𝑙 e − j 𝜔 𝑝 𝜏 𝑙 a R ( 𝜃 R ,𝑙 ) a T T ( 𝜃 T ,𝑙 ) x [ 𝑝, 𝑏 ] , (12) a T ( 𝜃 T ,𝑙 ) = (cid:104) 𝑒 j 𝜔 𝑐 𝜏 T , ( 𝜃 T ,𝑙 ) , . . . , e j 𝜔 𝑐 𝜏 T ,𝑁 T ( 𝜃 T ,𝑙 ) (cid:105) T ∈ C 𝑁 T (13)is the Tx array steering vector, with the Rx steering vector a R ( 𝜃 R ,𝑙 ) defined accordingly, 𝜔 𝑝 = 𝜋 𝑝 Δ 𝑓 , 𝜔 c = 𝜋 𝑓 c , ℎ 𝑙 is the gain of the 𝑙 -th path and η [ 𝑝, 𝑏 ] ∼ N C ( , 𝜎 𝜂 I 𝑁 R ) is the additivewhite Gaussian noise (AWGN). We write the signal model (11) as Y 𝑏 = ∑︁ 𝐿 − 𝑙 = ℎ 𝑙 C 𝑏 ( 𝜏 𝑙 , 𝜃 T ,𝑙 , 𝜃 R ,𝑙 ) + N 𝑏 , (14)where C 𝑏 ( 𝜏 𝑙 , 𝜃 T ,𝑙 , 𝜃 R ,𝑙 ) = a R ( 𝜃 R ,𝑙 ) a T T ( 𝜃 T ,𝑙 ) X 𝑏 diag ( a 𝜏 ( 𝜏 𝑙 )) ∈ C 𝑁 R × 𝑁 P (15) a 𝜏 ( 𝜏 ) = [ e − j 𝜔 𝑝 𝜏 , . . . , e − j 𝜔 𝑝𝑁 P 𝜏 ] T ∈ C 𝑁 P (16) Y 𝑏 = [ y [ 𝑝 , 𝑏 ] , . . . , y [ 𝑝 𝑁 P , 𝑏 ]] ∈ C 𝑁 R × 𝑁 P (17) X 𝑏 = [ x [ 𝑝 , 𝑏 ] , . . . , x [ 𝑝 𝑁 P , 𝑏 ]] ∈ C 𝑁 T × 𝑁 P (18) N 𝑏 = [ η [ 𝑝 , 𝑏 ] , . . . , η [ 𝑝 𝑁 P , 𝑏 ]] ∈ C 𝑁 R × 𝑁 P (19)Stacking the observations over 𝑁 B OFDM symbols we get Y = 𝐿 − ∑︁ 𝑙 = ℎ 𝑙 C ( 𝜏 𝑙 , 𝜃 T ,𝑙 , 𝜃 R ,𝑙 ) + N (20) where Y = [ Y T0 , . . . , Y T 𝑁 B − ] T (21) C ( 𝜏, 𝜃 T , 𝜃 R ) = [ C T0 ( 𝜏, 𝜃 T , 𝜃 R ) , . . . , C T 𝑁 B − ( 𝜏, 𝜃 T , 𝜃 R )] T (22) N = [ N T0 , . . . , N T 𝑁 B − ] . (23)Through (4), (7)-(8) and (20), we can see that the observations Y depend on the positionparameter vector ν , defined as ν = [ p T R , 𝛼 R , 𝜖 clk , h T0 , p T s , , h T1 , . . . , p T s ,𝐿 − , h T 𝐿 − ] T ∈ R 𝐿 + , (24)with h 𝑙 = [| ℎ 𝑙 | , arg ( ℎ 𝑙 )] T . C. Assumptions1) Reference signal structure:
In this work we consider the case where Tx uses a fixed beamcodebook f 𝑘 , 𝑘 = , . . . , 𝑀 T . This does not only simplify the optimization task, but also mightbe a practical limitation in a 5G system, with devices using a predefined set of beams fortransmission or reception.We also assume that the resource allocation R 𝑘 among the codebook beams and the powerallocation 𝛾 𝑘 [ 𝑝, 𝑏 ] among assigned REs, are fixed and therefore, optimizing R 𝑘 is not in thescope of our reference signal optimization task. The problem of designing a waveform the hasbeen addressed in [42]–[44], where the CRLB is optimized with respect to the resource allocationand additional constraints may be considered to avoid a unbalanced use of the spectrum.
2) Prior knowledge at Rx and Tx:
In many cases the Tx might have prior knowledge on ν , based on prior estimation in the reverse link, map information and known geographicaldistribution of the users. The prior information is encoded by the joint probability density function(pdf) 𝑝 ν ( ν ) . In the following, we examine how the Tx can expoit the prior information, so asto improve the ability to localize the Rx.The Rx, which aims to compute its position and orientation from the received signal, only hasknowledge on the clock offset’s distribution 𝑝 𝜖 clk , which we assume to be zero-mean Gaussianwith variance 𝜎 clk . We note that 𝜎 clk = and 𝜎 clk → ∞ correspond to perfect synchronizationand asynchronous operation, respectively. III. P
OSITION E RROR B OUND
The achievable positioning accuracy of the Rx can be characterized in terms of the hybridCRLB. For a parameter vector ν containing both deterministic and random paramters, thecovariance matrix C of any unbiased estimator ˆ ν of ν satisfies [45], [46] C − J − ν (cid:23) , (25)where (cid:23) denotes positive semi-definiteness and J ν ∈ R ( 𝐿 + )×( 𝐿 + ) is the hybrid Fisherinformation matrix (FIM) of ν . J ν is defined as J ν = J ( p ) ν + J ( o ) ν , (26)where J ( p ) ν = E ν 𝑟 [− 𝐷 ν ln 𝑝 ( ν 𝑟 )] , (27)accounts for the prior information and J ( o ) ν = E Y , ν 𝑟 [− 𝐷 ν ln 𝑝 ( Y | ν )] (28)accounts for the observation-related information, with ν 𝑟 representing the random parameters in ν . As 𝜖 clk is the only parameter with prior information at the Rx, it is straightforward to findthat, based on (24), the only non-zero entry of J ( p ) ν is (cid:2) J ( p ) ν (cid:3) , = / 𝜎 clk . (29)Since ν is observed under AWGN, the ( 𝑖, 𝑗 ) -th entry of the J ( o ) ν is (cid:104) J ( o ) ν (cid:105) 𝑖, 𝑗 = 𝜎 𝜂 𝑁 B ∑︁ 𝑏 = ∑︁ 𝑝 ∈P (cid:60) (cid:40) 𝜕 m H 𝑏 [ 𝑝 ] 𝜕𝜈 𝑖 𝜕 m 𝑏 [ 𝑝 ] 𝜕𝜈 𝑗 (cid:41) . (30)Using (4), (12) and (30), we can see that J ( o ) ν is independent of the value of 𝜖 clk . The SPEB isdefined as SPEB = tr ( E T J − ν E ) , (31)where E = [ e , e ] and e 𝑖 is the 𝑖 -th column of the identity matrix of the appropriate size. Theposition error bound (PEB) is defined as its square root. IV. B
EAM P OWER A LLOCATION O PTIMIZATION
For the reference signal optimization, we make use of the assumption that with large bandwisthand number of antennas the paths are asymptotically orthogonal [9], [12]. We note that the SPEBis a function of ν (cid:48) = [ p T R , 𝛼 R , | ℎ | , p T s , , | ℎ | , . . . , p T s ,𝐿 − , | ℎ 𝐿 − | T ] T ∈ R 𝐿 + , (32)that is, it is independent of the values of arg ( ℎ 𝑙 ) , 𝑙 = , . . . , 𝐿 − , and 𝜖 clk . Also, due to the innerproduct of the derivatives in (30), we can observe (see (9), (10) and (12)) that J is independentof 𝛽 𝑘 [ 𝑝, 𝑏 ] . In the following, we write J ν = J ν ( q , ν (cid:48) ) , with q = [ 𝑞 , . . . , 𝑞 𝑀 T ] ∈ R 𝑀 T , tostress that J ν is the hybrid FIM of ν , whose value depends on q and ν (cid:48) . Similarly, we writeSPEB = SPEB ( q , ν (cid:48) ) .We study how the Tx can optimize the beam power allocation q using its prior knowledgeon ν (cid:48) so as to enable higher positioning accuracy at the Rx. We choose the expected SPEB(ESPEB) ESPEB = E ν (cid:48) [ SPEB ( q , ν (cid:48) )] (33)as the performance metric. The following proposed methods can be easily adapted for otherobjectives, such as max ν (cid:48) SPEB ( q , ν (cid:48) ) . A. Problem formulation
The optimization problem in hand reads as: min q E ν (cid:48) [ SPEB ( q , ν (cid:48) )] s.t. q (cid:60) , T q ≤ , (34)where (cid:60) denotes element-wise inequality. In order to solve (34), one can employ a cubature rule[47], [48] with positive weights to approximate the expectation integral with a sum: E ν (cid:48) [ SPEB ( q , ν (cid:48) )] ≈ ∑︁ 𝑁 (cid:48) ν 𝑗 = 𝑝 𝑗 SPEB ( q , ν (cid:48) 𝑗 ) , (35)where ν (cid:48) 𝑗 and 𝑝 𝑗 > , 𝑗 = , . . . , 𝑁 ν (cid:48) are the cubature points and their corresponding weights,with 𝑁 ν (cid:48) being the number of cubature points. 𝑁 ν (cid:48) is determined by the dimension of ν (cid:48) andthe degree 𝑟 of the cubature . The cubature points and their weights are determined by the pdfof ν (cid:48) and 𝑟 . Then, (34) becomes min q ∑︁ 𝑁 ν (cid:48) 𝑗 = 𝑝 𝑗 SPEB ( q , ν (cid:48) 𝑗 ) s.t. q (cid:60) , T q ≤ . (36) A cubature rule has degree 𝑟 if it is exact for a (multivariate) polynomial of degree 𝑟 . In a similar fashion to [36], using the epigraph form of (36), we can show that it is equivalentto the following SDP: min q , B ,..., B 𝑁 ν (cid:48) ∑︁ 𝑁 ν (cid:48) 𝑗 = 𝑝 𝑗 tr ( B 𝑗 ) s.t. B 𝑗 E T E J ( q , ν (cid:48) 𝑗 ) (cid:23) , 𝑗 = , . . . , 𝑁 ν (cid:48) q (cid:60) , T q ≤ , (37)where B 𝑗 ∈ R × , 𝑗 = , . . . , 𝑁 ν are auxiliary variables of the SDP and (cid:23) denotes positivesemidefiniteness. The positivity requirement on the cubature weights is imposed to ensureconvexity of the objective in (37).The optimal vector q obtained with (37) may indicate that a very low power should be allocatedin the direction of the LOS path, which may lead to a missed detection of the LOS path at theRx. This can be avoided by ensuring that the excitation on directions around the LOS path is atleast a fraction 𝑞 th of the excitation in any other direction. To this end, for a given confidencelevel 𝜅 , we define 𝜃 ( 𝜅 ) T ,𝑙, min and 𝜃 ( 𝜅 ) T ,𝑙, max as the minimum and maximum AODs corresponding tothe two-dimensional (2D) Rx locations ( 𝑙 = ) or scatterer/reflector locations ( 𝑙 = , . . . , 𝐿 − )in the 𝜅 -confidence ellipse of the respective marginal. With a uniform grid of 𝑁 𝜃 possible AODs 𝜃 T ,𝑙,𝑚 within the interval [ 𝜃 ( 𝜅 ) T ,𝑙, min , 𝜃 ( 𝜅 ) T ,𝑙, max ] 𝜃 ( 𝜅 ) T ,𝑙,𝑚 = 𝜃 ( 𝜅 ) T ,𝑙, min + 𝑚 − 𝑁 𝜃 − 𝜃 ( 𝜅 ) T ,𝑙, max , 𝑚 = , . . . , 𝑁 𝜃 , (38)we define the excitation matrix A 𝑙 ∈ R 𝑁 𝜃 × 𝑀 T for the 𝑙 -th path as [ A 𝑙 ] 𝑚,𝑘 = | a T T ( 𝜃 ( 𝜅 ) T ,𝑙,𝑚 ) f 𝑘 | . (39)Finally, the excitation vector for the possible AODs of the 𝑙 -th path is A 𝑙 q . Finally, the vectorwith the excitation of the possible AODs associated with the 𝑙 -th path is A 𝑙 q . We augment (37)with the following linear constraints: A q (cid:60) 𝑞 th (cid:107) Aq (cid:107) ∞ 𝑁 𝜃 , (40)where A = [ A T0 , . . . , A T 𝐿 − ] T . We note that the constraints (40) can be equivalently expressedas A q (cid:60) 𝑞 th 𝑒 max 𝑁 𝜃 , Aq (cid:52) 𝑒 max 𝐿𝑁 𝜃 , (41)with 𝑒 max being an auxiliary optimization variable. The main challenge with the approach described above is that 𝑝 ν is a multidimensionalpdf. The number of auxiliary matrices B 𝑗 and corresponding positive semidefiniteness (PSD)constraints in (37) is equal to the number of cubature points. For known cubature rules [47], thenumber of points is lower bounded by ( 𝐿 + ) ( 𝑟 − )/ , which could result in very high complexityfor our optimization task, as the integrand is highly non-linear and a rule with 𝑟 ≥ is requiredfor an accurate approximation. B. Low-complexity sub-optimal solution1) Dimensionality reduction:
A way to circumvent the dimensionality challenge is to use asurrogate function which involves the expectation over a smaller set of parameters. To this end,we first note that e T 𝑖 J − e 𝑖 , 𝑖 = , , is a convex function of J and so is the SPEB as a sum ofconvex functions. Splitting ν (cid:48) into any couple of vectors ν and ν , we can write E ν [ SPEB ( q , ν )] = E ν (cid:2) tr ( E T J − ( q , ν (cid:48) ) E ) (cid:3) = E ν (cid:2) E ν | ν (cid:2) tr ( E T J − ( q , ν , ν ) E ) (cid:3) (cid:3) ( 𝑎 ) ≥ E ν (cid:2) tr ( E T ( E ν | ν [ J ( q , ν , ν )]) − E ) (cid:3) (42)where (a) follows from Jensen’s inequality. We choose ν = [ p T R , p T s , , . . . , p T s ,𝐿 − ] T , as theposition parameters are the ones determining the AODs, which in turn determine which beamsare relevant or not. One could optimize the lower bound on the ESPEB in (42), as described in(34)-(37), but the number of cubature points 𝑁 ν (cid:48) is still lower bounded by ( 𝐿 ) ( 𝑟 − )/ .
2) Power allocation as a weighted sum of per-path power allocation vectors:
Our aim isto reduce the complexity of the optimization problem in hand. We accomplish this by takingthe following heuristic approach: we compute a power allocation vector q 𝑙 , 𝑙 = , . . . , 𝐿 − ,considering the uncertainty regarding each path separately and then weight the resulting powerallocation vectors in order to minimize a lower bound on the ESPEB.More specifically, for the power allocation vector q we consider only the LOS path andneglect the NLOS paths and solve q = argmin q E p R (cid:2) tr ( E T ( E | ℎ | ,𝛼 R | p R [ J ν LOS ( q , p R , 𝛼 R , | ℎ |)]) − E ) (cid:3) s.t. A q (cid:60) 𝑞 th , LOS (cid:107) A q (cid:107) ∞ 𝑁 𝜃 , q (cid:60) , T q ≤ , (43)where J ν LOS represents the FIM for the parameter vector ν LOS = [ p T R , 𝛼 R , 𝜖 clk , h T0 ] T . Similarly to(40), the first constraint in (43) limits the ratio of power spent among possible LOS directions,with 𝑞 th , LOS being the corresponding minimum ratio. For the gain of the LOS path it is natural that 𝑝 ( h | p R ) = 𝑝 ( h | 𝑑 ) , with 𝑑 = (cid:107) p R (cid:107) , i.e. the distribution of the gain depends only on theTx-Rx distance, and the integration over the radial component 𝑑 and the angular component 𝜃 T , of p R can be carried out separately. Then, as shown in the Appendix, we can reformulate (43) asan SDP using a one-dimensional (1D) quadrature rule for the approximation of the expectationintegral over 𝜃 T , .For the power allocation vector q 𝑙 we consider only the 𝑙 -th NLOS path and assume thatthe Rx position and orientation are known and equal to their mean values ¯ p R and ¯ 𝛼 R . Thisis basically a bistatic radar setup, where the goal is the estimation of the point of incidence.Therefore, we obtain q 𝑙 by solving q 𝑙 = argmin q E p s ,𝑙 (cid:2) tr ( E T ( E | ℎ 𝑙 || p s ,𝑙 [ J NLOS ,𝑙 ( q , p s ,𝑙 , | ℎ 𝑙 |)]) − E ) (cid:3) s.t. q (cid:60) , T q ≤ , (44)where J NLOS ,𝑙 represent the FIM for the parameter vector ν NLOS ,𝑙 = [ p T s ,𝑙 , 𝜖 clk , h T 𝑙 ] T . Problem (44)can be solved employing a 2D cubature on p s ,𝑙 .Finally, we compute the optimal weights w ∈ R 𝐿 of q 𝑙 , 𝑙 = , . . . , 𝐿 − , by minimizing anapproximate lower bound on the ESPEB, obtained similarly to (42): w = argmin w (cid:48) E p R [ tr ( E T J − ( Qw (cid:48) , ¯ ν ) E )] s.t. A Qw (cid:48) (cid:60) 𝑞 th (cid:107) AQw (cid:48) (cid:107) ∞ 𝑁 𝜃 Qw (cid:48) (cid:60) , T Qw (cid:48) ≤ , (45)where, in order to further reduce the computational load, we have replaced E ν | p R [ J ( Qw (cid:48) , ν )] with its approximation J ( Qw (cid:48) , ¯ ν ) , with ¯ ν = E ν | p R [ ν ] and Q = [ q , . . . , q 𝐿 − ] . Finally, thebeam power allocation vector is q = Qw .V. C HANNEL AND P OSITION E STIMATION
In this section we present a novel algorithm for Rx position, orientation and clock offsetestimation. In the first step of the algorithm a gridless parameter estimation algorithm based on[41] is employed to recover the number paths and their respective TOAs, AODs and AOAs. Inthe second step, the recovered channel parameters are mapped to the position parameter vector ν . A. Channel parameter estimation
For our positioning purposes, we are not merely interested in denoising Y , but we wouldlike to recover the number of paths, along with their respective gains, TOAs, AODs and AOAs.Hence, we aim to solve the following optimization problem: min 𝐿 (cid:48) , { 𝜏 𝑙 ,𝜃 T ,𝑙 ,𝜃 R ,𝑙 ,ℎ 𝑙 } 𝐿 (cid:48)− 𝑙 = Λ ( R ) + 𝜒 (cid:107) h (cid:107) (46)where Λ ( R ) = (cid:107) R (cid:107) F (47)is the loss function, R = Y − ∑︁ 𝐿 (cid:48) − 𝑙 = ℎ 𝑙 C ( 𝜏 𝑙 , 𝜃 T ,𝑙 𝜃 R ,𝑙 ) (48)is the residual, 𝜒 is a regularization parameter and h = [ ℎ , . . . , ℎ 𝐿 (cid:48) − ] T . The penalty term (cid:107) h (cid:107) is included to make the channel representation more parsimonious; otherwise the numberof detected paths could grow arbitrarily so as to minimize the objective. As usual in sparserecovery setups, instead of a non-convex L0 norm penalty term, we use the L1 norm. We solveproblem (46) using the algorithmic framework of [41], termed as Alternating Descent ConditionalGradient Method (ADCGM), which is described in Alg. 1. We note that, for notational brevity,in (46)-(48) and in the following, we write R instead of R ( 𝐿 (cid:48) , { 𝜏 𝑙 , 𝜃 T ,𝑙 , 𝜃 R ,𝑙 } 𝐿 (cid:48) − 𝑙 = ) . Also, theresidual at iteration 𝑖 is denoted as R 𝑖 and the TOAs of the detected paths are stacked in thevector τ ( 𝑖 ) = [ 𝜏 ( 𝑖 ) , . . . , 𝜏 ( 𝑖 ) 𝐿 ( 𝑖 ) − ] ∈ R 𝐿 ( 𝑖 ) , where 𝐿 ( 𝑖 ) is the number of detected paths at iteration 𝑖 .The parameter vectors θ ( 𝑖 ) T and θ ( 𝑖 ) R are defined accordingly. The maximum number of iterationsis 𝐿 max and at each iteration a new path can be detected (Step 2) or previously detected pathscan be dropped (Step 4(b)). In the following, we describe steps 2 and 4 in detail.
1) Detection of a new potential path (Step 2):
In order to get the next potential path we haveto solve (49), which is non-convex and can be solved by discretizing the three-dimensional (3D)parameter space [ , 𝑇 CP ] × [− 𝜋, 𝜋 ) × [− 𝜋, 𝜋 ) to get an 𝑁 𝜏 × 𝑁 𝜃 T × 𝑁 𝜃 R -dimensional grid G .After computing the new potential source, we compare the correspoding objective with apredefined threshold 𝜁 > , which is a function of the noise variance 𝜎 𝜂 , the reference signal X and the desired false alarm probability 𝑃 fa . Algorithm 1
Channel parameter estimation with ADCGM input: { X 𝑏 } 𝑁 B 𝑏 = , Y , 𝜎 𝜂 , 𝑃 fa initialize: τ ( ) , θ ( ) T , θ ( ) R , h ( ) = [ ] do
1. Compute residual R 𝑖
2. Detect next potential path: 𝜏 ( 𝑖 ) , 𝜃 ( 𝑖 ) T , 𝜃 ( 𝑖 ) R = argmax ( 𝜏,𝜃 T ,𝜃 R )∈G (cid:12)(cid:12) tr ( R H 𝑖 C ( 𝜏, 𝜃 T , 𝜃 R )) (cid:12)(cid:12) (49)3. Update support: τ ( 𝑖 + ) = [( τ ( 𝑖 ) ) T , 𝜏 ( 𝑖 ) ] , θ ( 𝑖 + ) T = [( θ ( 𝑖 ) T ) T , 𝜃 ( 𝑖 ) T ] , θ ( 𝑖 + ) R = [( θ ( 𝑖 ) R ) T , 𝜃 ( 𝑖 ) R ]
4. Coordinate descent on non-convex objective: for 𝑖 = to 𝑁 cd do (a) Compute gains: h ( 𝑖 + ) = argmin h Λ ( R ) + 𝜒 (cid:107) h (cid:107) (50)(b) Prune support: { τ , θ T , θ R , h } ( 𝑖 + ) = prune ({ τ , θ T , θ R , h } ( 𝑖 + ) ) (c) Locally improve support: { τ , θ T , θ R } ( 𝑖 + ) = local _ descent ({ τ , θ T , θ R , h } ( 𝑖 + ) ) end for 𝑖 = 𝑖 + while 𝑘 < 𝐿 max and (cid:12)(cid:12) tr ( R H 𝑖 C ( 𝜏 ( 𝑖 ) , 𝜃 ( 𝑖 ) T , 𝜃 ( 𝑖 ) R ) (cid:12)(cid:12) > 𝜁
2) Coordinate descent (Step 4):
In this algorithmic step we iteratively perform 3 sub-stepsfor a fixed number of 𝑁 cd iterations:(a) We update the gains solving (50), keeping the other path parameters fixed. The regularizationparameter 𝜒 determines the accuracy-sparsity trade-off.(b) We prune the paths whose gain is effectively zero: the 𝑙 -th path is pruned if | ℎ 𝑙 | / 𝜁 < max 𝑙 = ,...,𝐿 ( 𝑖 ) − | ℎ 𝑙 | , where < 𝜁 (cid:28) .(c) For the local descent step we perform truncated Newton steps for each path and each parameter sequentially: 𝜏 ( 𝑖 + ) 𝑙 ← 𝜏 ( 𝑖 + ) 𝑙 − sgn ( 𝜕 Λ / 𝜕𝜏 ( 𝑖 + ) 𝑙 ) 𝑠 ( 𝑖 + ) 𝜏,𝑙 (51) 𝜃 ( 𝑖 + ) T ,𝑙 ← 𝜃 ( 𝑖 + ) T ,𝑙 − sgn ( 𝜕 Λ / 𝜕𝜃 ( 𝑖 + ) T ,𝑙 ) 𝑠 ( 𝑖 + ) 𝜃 T ,𝑙 (52) 𝜃 ( 𝑖 + ) R ,𝑙 ← 𝜃 ( 𝑖 + ) R ,𝑙 − sgn ( 𝜕 Λ / 𝜕𝜃 ( 𝑖 + ) R ,𝑙 ) 𝑠 ( 𝑖 + ) 𝜃 R ,𝑙 (53)where 𝑠 ( 𝑖 + ) 𝜏,𝑙 = min (cid:18)(cid:12)(cid:12)(cid:12)(cid:16) 𝜕 Λ /( 𝜕𝜏 ( 𝑖 + ) 𝑙 ) (cid:17) − 𝜕 Λ / 𝜕𝜏 ( 𝑖 + ) 𝑙 (cid:12)(cid:12)(cid:12) , 𝑁 CP 𝑇 s ( 𝑁 𝜏 − ) (cid:19) 𝑠 ( 𝑖 + ) 𝜃 T ,𝑙 = min (cid:18)(cid:12)(cid:12)(cid:12)(cid:16) 𝜕 Λ /( 𝜕𝜃 ( 𝑖 + ) T ,𝑙 ) (cid:17) − 𝜕 Λ / 𝜕𝜃 ( 𝑖 + ) T ,𝑙 (cid:12)(cid:12)(cid:12) , 𝜋𝑁 𝜃 T − (cid:19) 𝑠 ( 𝑖 + ) 𝜃 R ,𝑙 = min (cid:18)(cid:12)(cid:12)(cid:12)(cid:16) 𝜕 Λ /( 𝜕𝜃 ( 𝑖 + ) R ,𝑙 ) (cid:17) − 𝜕 Λ / 𝜕𝜃 ( 𝑖 + ) R ,𝑙 (cid:12)(cid:12)(cid:12) , 𝜋𝑁 𝜃 R − (cid:19) are the step sizes, with 𝑇 s = 𝑁 Δ 𝑓 . We note that we limit the maximum step size for eachof the parameters to be equal to half of the corresponding grid bin size, in order to avoidconvergence problems near inflection points of the loss function. B. Mapping to position parameters
Having an estimate ˆ˜ ν of the channel parameter vector ˜ ν defined as ˜ ν = [ 𝜏 , 𝜃 T , , 𝜃 R , , . . . , 𝜏 ˆ 𝐿 − , 𝜃 T , ˆ 𝐿 − , 𝜃 R , ˆ 𝐿 − ] T , (54)where ˆ 𝐿 is the estimated number of paths, and choosing the strongest path as the LOS path,we estimate the position parameter vector ν employing the EXIP as in [13], with a slightmodification to include the prior information on the clock offset. To this end, we intend to solve argmin ν ( ˆ˜ ν − 𝑓 ( ν )) T J ˆ˜ ν ( ˆ˜ ν − 𝑓 ( ν )) + ( 𝜖 clk / 𝜎 clk ) , (55)where J ˆ˜ ν is the channel parameter FIM and 𝑓 : R 𝐿 + → R 𝐿 is the mapping from position tochannel parameters, determined by (4), (7)-(8).We note that false alarms, that is falsely detected paths, can have severe impact on positionestimation. Therefore, we apply the following two criteria to filter them out: • A single-bounce NLOS path and a LOS path always form a triangle, as can be seen inFig. 1. Therefore, for the formation of a triangle to be possible, a single-bounce NLOSpath must satisfy Δ 𝜃 T ,𝑙 · Δ 𝜃 R ,𝑙 < , 𝑙 = , . . . , ˆ 𝐿 − , (56) where Δ 𝜃 T ,𝑙 = 𝜃 T ,𝑙 − 𝜃 T , and Δ 𝜃 R ,𝑙 = 𝜃 R ,𝑙 − 𝜃 R , , with Δ 𝜃 T ,𝑙 and Δ 𝜃 R ,𝑙 ∈ [− 𝜋, 𝜋 ) . Thereforeif the 𝑙 -th path, 𝑙 = , . . . , ˆ 𝐿 − , does not satisfy (56), it is dropped. • Combined with the LOS path, each NLOS path can provide an estimate of 𝜖 clk : 𝜖 clk ,𝑙 = 𝜏 𝑙 sin ( Δ 𝜃 R ,𝑙 − Δ 𝜃 T ,𝑙 ) − 𝜏 ( sin ( Δ 𝜃 R ,𝑙 ) − sin ( Δ 𝜃 T ,𝑙 )) sin ( Δ 𝜃 R ,𝑙 − Δ 𝜃 T ,𝑙 ) − ( sin ( Δ 𝜃 R ,𝑙 ) − sin ( Δ 𝜃 T ,𝑙 )) , (57)With 𝜁 ,𝑎 > and 𝜁 ,𝑏 > being predefined probability thresholds for 𝜖 clk values, if 𝑝 ( 𝜖 clk ,𝑙 ) < 𝜁 ,𝑎 or 𝑝 ( 𝜖 clk ,𝑙 ) < 𝜁 ,𝑏 𝑝 clk , max , the path is filtered out, with 𝑝 clk , max = max 𝑙 = ,..., ˆ 𝐿 − 𝑝 ( 𝜖 clk ,𝑙 ) .Replacing ˆ˜ ν with ˆ˜ ν (cid:48) , which contains only the remaining paths, we solve (55) with theLevenberg-Marquardt algorthm [49], [50]. For the initial point ν ( ) we compute 𝜖 ( ) clk = (cid:205) 𝑙 | ℎ 𝑙 | 𝜖 clk ,𝑙 (cid:205) 𝑙 | ℎ 𝑙 | (58) p ( ) R = 𝑐 ( 𝜏 − 𝜖 ( ) clk ) u ( 𝜃 T , ) (59) 𝛼 ( ) R = 𝜃 T , + 𝜋 − 𝜃 R , (60) p ( ) s ,𝑙 = tan ( 𝜃 R ,𝑙 + 𝛼 ( ) R ) 𝑝 ( ) R ,𝑥 − 𝑝 ( ) R ,𝑦 tan ( 𝜃 R ,𝑙 + 𝛼 ( ) R ) cos 𝜃 T ,𝑙 − sin 𝜃 T ,𝑙 u ( 𝜃 T ,𝑙 ) , 𝑙 = , . . . , ˆ 𝐿 (cid:48) , (61)where ˆ 𝐿 (cid:48) is the number of remaining estimated paths.VI. N UMERICAL R ESULTS
A. Simulation setup1) Geometric setup and prior information at the Tx:
For the evaluation of the power allocationand the position estimation algorithms we consider the setup shown in Fig. 2. The Tx is equippedwith a ULA with 𝑁 T = antennas. In order to be able to discriminate all possible AOAs, theRx has a UCA with 𝑁 R = antennas. With the Rx being equipped with a UCA, the SPEB isindependent of the orientation 𝛼 R .We consider NLOS paths resulting from single-bounce reflections. The phases of the complexpath gains are uniformly distributed over [− 𝜋, 𝜋 ) and their magnitudes are given by | ℎ 𝑙 | = 𝑐 /( 𝜋 𝑓 c (cid:107) p R (cid:107) ) , 𝑙 = , √ 𝜌 𝑙 𝑐 /( 𝜋 𝑓 c ( (cid:107) p s ,𝑙 (cid:107) + (cid:107) p R − p s ,𝑙 (cid:107) )) , 𝑙 ≠ , (62) − − p T p R p s , p s , p s , Fig. 2. Prior knowledge at the Tx for simulation results. where 𝜌 𝑙 is the reflection coefficient and 𝜆 c = 𝑐 / 𝑓 c . The prior knowledge at the Tx is describedby N ( µ , C ) , where µ = [ ¯ p T R , ¯ p T s , , ¯ 𝜌, ¯ p T s , , ¯ 𝜌, ¯ p T s , , ¯ 𝜌 ] T ∈ R (63) C = C , C , C , C , C T0 , C , 𝜎 𝜌 C T0 , C , 𝜎 𝜌 C T0 , C ,
00 0 𝜎 𝜌 ∈ R × (64)with ¯ p R = m , C , = /√ I m ¯ p s , = . m , C , = .
48 00 1 m , C , = .
45 00 0 m ¯ p s , = . − m , C , = .
34 00 1 m , C , = .
64 00 0 m ¯ p s , = . m , C , = . m , C , = . m ¯ 𝜌 = − , 𝜎 𝜌 = . Samples from this distribution are depicted in Fig. 2.
2) System Parameters:
For the waveform we set 𝑓 c =
38 GHz , 𝑁 = , 𝑁 B = , P = {− , . . . , − , . . . , } and Δ 𝑓 ( max (P) − min (P)) (≈ 𝐵 ) =
120 MHz . The resources are as-signed to the beams in an interleaved and staggered manner, i.e. R 𝑘 = {( 𝑘 + 𝑏 + 𝑖 𝑀 T , 𝑏 ) | 𝑖 ∈ Z , 𝑏 = , . . . , 𝑁 B : 𝑘 + 𝑏 + 𝑖 𝑀 T ∈ P } . The power of each beam is distributed uniformly amongits resources, i.e., 𝛾 𝑘 [ 𝑝, 𝑏 ] = /|R 𝑘 | . The noise variance is 𝜎 𝜂 = . ( 𝑛 Rx + 𝑁 ) 𝑁 Δ 𝑓 , where 𝑁 = −
174 dBm Hz − is the noise power spectral density per dimension and 𝑛 Rx = isthe Rx noise figure. The standard deviation of the clock offset is 𝜎 clk = /( 𝑁 Δ 𝑓 ) , so that 𝑐𝜎 clk ≈ .
88 m . We use a DFT beam codebook: f 𝑘 = (cid:2) , e − j 𝜋𝑁 T 𝑘 , . . . , e − j 𝜋𝑁 T ( 𝑁 T − ) 𝑘 (cid:3) , 𝑘 = , . . . , 𝑀 T = 𝑁 T . (65)
3) Benchmark for beam power allocation:
In order to fairly evaluate our power allocationstrategies, we set as benchmark the uniform power allocation to beams exciting useful directions.For a given confidence level 𝜅 we get a grid of AODs for each path as in (38) and compute theset of useful beams as B ( 𝜅 ) uni = ∪ 𝐿 − 𝑙 = ∪ 𝑁 𝜃 𝑚 = (cid:110) argmax 𝑘 = ,...,𝑁 T | a T T ( 𝜃 ( 𝜅 ) T ,𝑙,𝑚 ) f 𝑘 | (cid:111) . (66)The power allocation vector q is 𝑞 𝑘 = /|B ( 𝜅 ) uni | , 𝑘 ∈ B ( 𝜅 ) uni , 𝑘 ∉ B ( 𝜅 ) uni . (67) B. Power allocation strategies and position estimation algorithm parameters
The power allocation strategies and their corresponding parametrizations that we consider forour simulation results are as follows: • opt. unconstr. : Solution of (37). The number of points of known cubatures of 5th degree(in order to ensure a sufficiently dense sampling of the support of the distribution) withpositive weights is + · = , which incurs prohibitive computational complexity.Instead, we draw = random samples (as many as the lower bound for any cubature)from the joint -dimensional distribution. • opt. constr. : Solution of (37) with random samples from the joint -dimensionaldistribution and additional constraints (40), with 𝜅 = . , 𝑞 th = −
10 dB and 𝑁 𝜃 = . − − − L O S N L O S N L O S N L O S (a) opt. unconstr. − − − L O S N L O S N L O S N L O S (b) opt. constr. − − − L O S N L O S N L O S N L O S (c) opt. reduced − − − L O S N L O S N L O S N L O S (d) subopt. − − − L O S N L O S N L O S N L O S (e) uni 0.60 − − − L O S N L O S N L O S N L O S (f) uni 0.90Fig. 3. Beam patterns | a T T ( 𝜃 T ) f 𝑘 √ 𝑞 𝑘 | , 𝑘 = , . . . , 𝑀 T , for different power allocation strategies. • opt. reduced : Solution of the minimization of the lower bound on ESPEB (42) with = random samples from the joint -dimensional distribution and additional constraints (40),with 𝜅 = . , 𝑞 th = −
10 dB and 𝑁 𝜃 = . • subopt. : Solution of (43)-(45), with -point cubatures for the involved 2D marginals, 𝜅 = . , 𝑞 th , LOS = − , 𝑞 th = −
10 dB and 𝑁 𝜃 = . • uni 𝜅 : Uniform power allocation to useful directions, according to (66)-(67), with 𝜅 = { . , . } and 𝑁 𝜃 = . We note that choosing 𝜅 = . as for the other strategiesresults in performance degradation; hence, results for this value are not icluded.The beampatterns of the power allocation strategies for the considered prior knowledge are shownin Fig. 3. We observe in Figs. 3(a)-(d) that for the optimized power allocation strategies, mostof the available power is spent on beams illuminating NLOS paths. When 𝜎 clk is very small (i.e.when the synchronization error is very small), having only the delay measurement of the LOSsuffices to determine the distance between the base station (BS) and the user equipment (UE).However, as 𝜎 clk increases, neither the LOS nor the NLOS provide individually information about the BS-UE distance. In these cases, they are the differences between delays that are informative,and this implies that several paths (not only one) have to be illuminated with sufficient powerbecause if there is a large power unbalance between rays, then the delay differences will not beestimated precisely. Comparing Fig. 3(a) with Figs. 3(b)-(d), we see that when the constraints (40)are not applied, the power allocation to NLOS components is more significant, with the powerinvested to less likely LOS directions being very low. From Figs. 3(b) and (c), we can seethat the impact of the dimensionality reduction (42) is the reduction of the power spent on the2nd NLOS path. This is explained by the fact that the fading of the path gains is not takeninto account; hence, for the mean values of the path gains, more power is spent on the pathsthat offer more useful position information. Also, in Fig. 3(d) we observe that our suboptimalapproach allocates almost no power to the 2nd NLOS path, as in the last step where all paths areconsidered jointly, only the receiver’s location uncertainty and the mean scatterers’/reflectors’locations are taken into account; for this setup, the information offered by the 1st NLOS path ismore useful and therefore most of the available power is allocated for its illumination. For theuniform allocation, higher confidence values lead to activation of more beams and spreading ofthe avaiable power to more directions.Regarding the position estimation algorithm parameters, we set 𝑁 𝜏 = 𝑁 P , 𝑁 𝜃 T = 𝑁 T , 𝑁 𝜃 R = 𝑁 R , 𝑃 fa = . , 𝜁 is pre-trained for the given 𝑃 fa and power allocation strategy, 𝜁 = −
35 dB , 𝑁 cd = , 𝐿 max = , 𝜒 = 𝜎 𝜂 √︁ ( 𝑁 T + 𝑁 R ) |P | 𝑁 B 𝑃 RE / 𝑁 T (chosen according to [51]), 𝜁 ,𝑎 = − and 𝜁 ,𝑏 = − . C. Performance vs SNR for fixed geometry
We fix the geometry and the reflection coefficients to their mean value µ in (63) to examinethe performance of the position estimation algorithm as a function of the Tx power. For thepower allocation strategies described in Sec. VI-B, in Fig. 4, we plot the position root meansquare error (RMSE) E η ,𝜖 clk [(cid:107) ˆ p R − p R (cid:107) ] and PEB as functions of the average power per resourceelement 𝑃 RE = 𝑃 tot /( 𝑁 B 𝑁 P ) , with ˆ p R being the position estimate. We note that the average Txpower 𝑃 T is related to 𝑃 RE as 𝑃 T = 𝑃 RE 𝑁 P / 𝑁 .We can see that the bound is attained for all power allocation strategies. Regarding uniformpower allocation, the distance of the RMSE from the bound for low Tx power is attributedto the fact that, although the LOS path is detected, the probability of detection for the NLOSis small. With only the LOS path being detected, the clock offset cannot be resolved and the − − − − − − − − P RE in dBm P o s iti on R M S E i n m opt. unconstr.opt. constr.opt. reducedsubopt.uni 0.60uni 0.90 Fig. 4. Position RMSE (solid lines) and PEB (dashed lines) vs Tx power for different power allocation strategies. resulting position RMSE approaches the standard deviation of the clock offset 𝑐 · 𝜎 clk ≈ .
88 m .Among the two considered configurations ( 𝜅 = . and 𝜅 = . ), the former has slightly betterperformance, as the available power is more concentrated to the true location of the Rx and thereflectors. But, as we will see later on, this comes with a cost, when the uncertainty about thegeometry is considered.The optimized allocation strategies result in similar PEBs and offer significant improvementcompared to the uniform, with a gain of to for the same localization accuracy. The lowestPEB is attained by "opt. unconstr.", but the RMSE converges to the PEB for larger 𝑃 RE , comparedto the other strategies. The reason for this behavior is that, as can be observed in Fig. 3(a), only asmall fraction of power is spent in the LOS direction and the Tx power required for the LOS pathto be detected is larger. When the LOS path is missed, the first arriving NLOS path is treated asLOS by the algorithm, resulting in a large position error. Due to the constraints (40), the rest ofthe proposed strategies ("opt. constr.", "opt. reduced" and "subopt.") allocate more power to theLOS, enabling to attain the PEB at lower values of 𝑃 RE , with only a small performance penatly.The RMSE of "opt. reduced" converges slightly faster to the bound compared to "opt. constr.",as slightly more power is allocated to the LOS path. The "subopt." allocation exhibits the mostrobust performance, as the LOS path can be detected for much lower Tx power values. . . . . . . . . . k ˆ p R − p R k in m E m p i r i ca l C D F opt. unconstr.opt.constr.opt. reducedsubopt.uni 0.60uni 0.90 Fig. 5. Empirical cdf of (cid:107) ˆ p R − p R (cid:107) for different power allocation strategies.TABLE IP ERCENTILES OF THE CDF OF THE POSITION ERROR IN m FOR DIFFERENT POWER ALLOCATION STRATEGIES .50% 90% 95% 99%opt. unconstr. 0.22 0.93 29.55 72.25opt. constr. 0.21 0.59 0.78 1.32opt. reduced 0.21 0.57 0.76 1.25subopt. 0.21 0.65 0.84 1.45uni 0.60 0.31 0.91 1.30 20.83uni 0.90 0.30 0.86 1.10 1.96
D. Performance with random geometry
The results in Fig. 4 and the corresponding discussion were useful in examining the behaviorof the position estimation algorithm, but do not provide a complete characterization of theperformance of the power allocation strategies. To better evaluate their performance, for 𝑃 RE = and the rest of the system parameters as described in Sec. VI-A2, we plot in Fig. 5the cumulative distribution function (cdf) of the position error (cid:107) ˆ p R − p R (cid:107) , which is computeddrawing samples from (63)-(64). A summary of the percentiles of the distribution of the positionerror is provided in Table I.We can observe in Fig. 5 and Table I that "opt. reduced" and "opt. constr." achieve the bestperformance. The latter is slightly worse at higher percentiles, as more points would be requiredfor a more accurate approximation of the expectation in the corresponding optimization problem. In spite of the lower computations cost of the "subopt." allocation, its performance degradationis almost unnoticeable. On the other hand, the "opt. unconstr." approach, although attainingalmost the same median error as the other optimized strategies, has much lower accuracy forhigher percentiles. This is attributed to the low power spent in the direction around the LOSpath, resulting in low probability of detection of the LOS. Compared to the best of the uniformallocations, the "opt. reduced" power allocation offers a position error reduction of 30%, 34%,31% and 36% at the 50%, 90%, 95% and 99% percentile, respectively.Regarding the uniform allocations, we can see that spreading the power to a reduced set ofbeams ("uni 0.60") might result in better positioning accuracy for some geometry realizations,as seen for example in Fig. 4, but it significantly deteriorates the performance for other possiblerealizations. This explains the higher values of position errors at the upper percentiles of thecorresponding cdf.
E. Power allocation as a function of 𝜎 clk We now examine the effect of 𝜎 clk on the power allocation. First, similar to (66), we definethe set of LOS-illuminating beams as B ( 𝜅 ) LOS = ∪ 𝑁 𝜃 𝑚 = (cid:110) argmax 𝑘 = ,...,𝑁 T | a T T ( 𝜃 ( 𝜅 ) T , ,𝑚 ) f 𝑘 | (cid:111) . (68)and the fraction of power spent on them as 𝑞 LOS = ∑︁ 𝑘 ∈B ( 𝜅 ) LOS 𝑞 𝑘 . (69)In Fig. 6(a) we plot 𝑞 LOS as a function of 𝜎 clk for the power allocation strategies "opt. unconstr.","opt. constr.", "subopt" and "uni 0.90", for 𝑁 R = { , } , 𝑃 RE = , 𝜅 = . and the restof the system parameters as described in Sec. VI-A2; in Fig. 6(b) we plot the corresponding E [ PEB ] . We can see in Fig. 6(a) that for very low values of 𝜎 clk , equivalent to almost perfectTx-Rx synchronization, it is optimal to spend almost all the available power on LOS-illuminatingbeams. As 𝜎 clk increases, 𝑞 LOS decreases rapidly for both optimized allocation strategies, untilit saturates at a relatively low value. This is explained as follows: The clock offset decreases theamount of range information provided by the LOS path and the larger standard deviation of theclock offset, the more significant the decrease. Hence, as 𝜎 clk increases, the ranging informationprovided by the NLOS paths becomes more significant and, therefore, more power is spenton them. Nevertheless, the saturation occurs because the measurement of the LOS AOD offers . . . . . . . c σ clk in m q L O S opt. unconstr., N R = opt. unconstr., N R = opt. constr., N R = opt. constr., N R = subopt., N R = subopt., N R = uni 0.90, N R = uni 0.90, N R = (a) 𝑞 LOS vs. 𝜎 clk . . . . . c σ clk in m E [ P E B ] i n m (b) E [ PEB ] vs. 𝜎 clk Fig. 6. Fraction of power allocated to LOS-illuminating beams 𝑞 LOS and E [ PEB ] as functions of 𝜎 clk . significant information in the orthogonal direction, which is reduced when 𝑞 LOS is decreased. Thesaturation value for "opt. constr." is higher due to the additional constraints on LOS illumination.Also, we observe that the transition from high to low 𝑞 LOS values is slower for 𝑁 R = . This isattributed to the fact that NLOS paths offer rank-1 position information, whose intensity dependson the quality of the TOA, AOD and AOA measurements combined [12]. For 𝑁 R = , the qualityof the AOA measurement is lower; therefore, the intensity of the ranging information from theNLOS paths is smaller, compared to 𝑁 R = , and becomes significant for larger values of 𝜎 clk .In Fig. 6(b) we see that E [ PEB ] increases with increasing 𝜎 clk , until it saturates at a valuedependent on the power allocation strategy and the system configuration ( 𝑁 R = { , } ). As 𝜎 clk increases the reduction of ranging information from the LOS path cannot be complementedby ranging information from the NLOS paths (even with optimized power allocation), resultingin a larger error. In the saturation region the ranging information from the LOS path becomesnegligible compared to the clock offset-independent part of ranging information offered by thecombination of NLOS paths with the LOS path.VII. C ONCLUSION
Optimal power allocation on a beam codebook for single-anchor localization and lower-complexity suboptimal alternatives have been considered under imperfect Tx-Rx synchronization. A channel and position estimation method has also been proposed. Numerical results showthat our suboptimal power allocation approach offers a good balance between performance andcomplexity, as the significant complexity reduction for the computation of the power allocationincurs only a very small performance penalty. Our analysis has shown that even for low valuesof the clock offset standard deviation it is optimal according to the CRLB to allocate most of theavailable power to scatterer/reflector illuminating beams to recover necessary range information.We have also shown that guaranteeing a minimum amount of power spent on LOS-illuminatingbeams, can be beneficial when the actual position estimation is considered, as it ensures that LOSpath is detected with a high probability. The proposed position estimation algorithm attains thecorresponding CRLB for all considered power allocation strategies, benefiting from the gridlessparameter estimation, which avoids the appearance of spurious paths due to grid mismatch, whilealso filtering out noisy detected paths exploiting the information on the clock offset carried bysingle-bounce-NLOS paths. A
PPENDIX P OWER A LLOCATION FOR THE
LOS P
ATH
Here we show how to formulate (43) as an SDP using only a 1D quadrature rule for theapproximation of the expectation over 𝜃 T , . This is accomplished in two steps: • In the first step we show that the integration over 𝑑 and 𝜃 T , can be carried out separately; • in the second step, after averaging over 𝑑 , we exploit the form of the resulting functionof 𝜃 T , and formulate the problem as an SDP.We write E 𝑑 ,𝜃 T , [·] instead of E p R [·] . Also, for notational brevity we write ¯ J = E 𝛼 R , h | 𝑑 ,𝜃 T , [ J ν LOS ( q , 𝑑 , 𝜃 T , , 𝛼 R , h )] . (70)We index the elements of ¯ J with the pair of parameters to which they correspond.First, after some algebra we find that tr ( E T ¯ J − E ) = 𝑐 ¯ 𝐽 𝜏 ,𝜏 − ¯ 𝐽 𝜏 ,𝜃 T , ¯ 𝐽 𝜃 T , ,𝜃 T , + 𝑑 ¯ 𝐽 𝜃 T , ,𝜃 T , − ¯ 𝐽 𝜏 ,𝜃 T , ¯ 𝐽 𝜏 ,𝜏 + 𝑐 𝜎 clk , (71)where ¯ 𝐽 𝑎,𝑏 = E 𝛼 R , h | 𝑑 ,𝜃 T , [ 𝐽 𝑎,𝑏 ] (72) 𝐽 𝑎,𝑏 = 𝜎 𝜂 𝑁 B ∑︁ 𝑏 = ∑︁ 𝑝 ∈P (cid:60) (cid:26) 𝜕 m H 𝑏 [ 𝑝 ] 𝜕𝑎 𝜕 m 𝑏 [ 𝑝 ] 𝜕𝑏 (cid:27) , (73) with 𝑎, 𝑏 ∈ { 𝑑 , 𝜃 T , } . We can show that 𝐽 𝑎,𝑏 , 𝑎, 𝑏 ∈ { 𝑑 , 𝜃 T , } , are independent of 𝛼 R and thephase of ℎ . Hence, they can be expressed as ¯ 𝐽 𝑎,𝑏 = E h | 𝑑 ,𝜃 T , [ 𝐽 𝑎,𝑏 ( q , 𝜃 T , , | ℎ ( 𝑑 ) | )] = E h | 𝑑 ,𝜃 T , [| ℎ ( 𝑑 ) | 𝑗 𝑎,𝑏 ( q , 𝜃 T , )] = 𝑔 ( 𝑑 ) 𝑗 𝑎,𝑏 ( q , 𝜃 T , ) , (74)where 𝑔 ( 𝑑 ) = E h | 𝑑 [| ℎ ( 𝑑 ) | ] and 𝑗 𝑎,𝑏 ( q , 𝜃 T , ) = 𝐽 𝑎,𝑏 ( q , 𝜃 T , , | ℎ ( 𝑑 ) | )/| ℎ ( 𝑑 ) | is a func-tion of q and 𝜃 T , . For the second equality in (74), we used the fact that 𝐽 𝑎,𝑏 can be expressedas the product of two terms, one dependent on the gain magnitude and the other on q and 𝜃 T , .We can then rewrite (71) as tr ( E T ¯ J − E ) = 𝑔 ( 𝑑 ) (cid:18) 𝑐 𝐼 𝜏 ( q , 𝜃 T , ) + 𝑑 𝐼 𝜃 T , ( q , 𝜃 T , ) (cid:19) + 𝑐 𝜎 clk , (75)where 𝐼 𝜏 ( q , 𝜃 T , ) = 𝑗 𝜏 ,𝜏 ( q , 𝜃 T , ) − 𝑗 𝜏 ,𝜃 T , ( q , 𝜃 T , ) 𝑗 𝜃 T , ,𝜃 T , ( q , 𝜃 T , ) (76) 𝐼 𝜃 T , ( q , 𝜃 T , ) = 𝑗 𝜃 T , ,𝜃 T , ( q , 𝜃 T , ) − 𝑗 𝜏 ,𝜃 T , ( q , 𝜃 T , ) 𝑗 𝜏 ,𝜏 ( q , 𝜃 T , ) . (77)It is then apparent from the form of the function in (75) that integration over 𝑑 and 𝜃 T , canbe carried out separately.For the second step, taking the expectation over 𝑑 and defining ¯ 𝑔 ( 𝜃 T , ) = / E 𝑑 | 𝜃 T , [ / 𝑔 ( 𝑑 )] (78) ¯ 𝑑 ( 𝜃 T , ) = √︄ E 𝑑 | 𝜃 T , (cid:20) ¯ 𝑔 ( 𝜃 T , ) 𝑔 ( 𝑑 ) 𝑑 (cid:21) (79)we get E 𝑑 | 𝜃 T , [ tr ( E T ¯ J − E )] = 𝑔 ( 𝜃 T , ) (cid:18) 𝑐 𝐼 𝜏 ( q , 𝜃 T , ) + (cid:0) ¯ 𝑑 ( 𝜃 T , ) (cid:1) 𝐼 𝜃 T , ( q , 𝜃 T , ) (cid:19) + 𝑐 𝜎 clk . (80)Comparing (80) to (71), we can conclude that, in order to be able to formulate the problem ina convex form, E 𝑑 | 𝜃 T , [ tr ( E T ¯ J − E )] can be expressed as E 𝑑 | 𝜃 T , [ tr ( E T ˇ J − E )] = tr ( E T J − ν LOS ( q , ¯ 𝑑 ( 𝜃 T , ) , 𝜃 T , , ˇ 𝛼 R , √︁ ¯ 𝑔 ( 𝜃 T , ) e j 𝛽 𝑔 ) E ) , (81)where ˇ 𝛼 R and 𝛽 𝑔 can be chosen arbitrarily, since they do not have an impact on the objective.Finally, using (81) and the identity E 𝑑 ,𝜃 T , (cid:2) tr ( E T ¯ J E ) (cid:3) = E 𝜃 T , [ E 𝑑 | 𝜃 T , [ tr ( E T ¯ J − E )]] , (82) we can employ a 1D quadrature rule to approximate the expectation integral over 𝜃 T , to get thefollowing SDP: min q , B ,..., B 𝑁𝜃 T , ∑︁ 𝑁 𝜃 T , 𝑗 = 𝑝 𝑗 tr ( B 𝑗 ) s.t. q ≥ , T q ≤ , B 𝑗 E T E J ν LOS ( q , ¯ 𝑑 ( 𝜃 T , , 𝑗 ) , 𝜃 T , , 𝑗 , ˇ 𝛼 R , √︁ ¯ 𝑔 ( 𝜃 T , , 𝑗 ) e j 𝛽 𝑔 ) (cid:23) , 𝑗 = , . . . , 𝑁 𝜃 T , . (83)R EFERENCES [1] S. Mumtaz, J. Rodriguez, and L. Dai, mmWave Massive MIMO: A Paradigm for 5G . Academic Press, 2017.[2] A. L. Swindlehurst, E. Ayanoglu, P. Heydari, and F. Capolino, “Millimeter-wave massive MIMO: the next wirelessrevolution?”
IEEE Commun. Mag. , vol. 52, no. 9, pp. 56–62, Sep. 2014.[3] F. Wen, H. Wymeersch, B. Peng, W. P. Tay, H. C. So, and D. Yang, “A survey on 5G massive MIMO localization,”
DigitalSignal Processing , vol. 94, pp. 21 – 28, Nov. 2019, Special Issue on Source Localization in Massive MIMO.[4] J. A. del Peral-Rosado, R. Raulefs, J. A. López-Salcedo, and G. Seco-Granados, “Survey of cellular mobile radio localizationmethods: From 1G to 5G,”
IEEE Commun. Surveys Tuts. , vol. 20, no. 2, pp. 1124–1148, 2nd quarter 2018.[5] R. Keating, M. Säily, J. Hulkkonen, and J. Karjalainen, “Overview of positioning in 5G new radio,” in
Proc. 16th Int.Symposium on Wireless Commun. Systems (ISWCS) , Aug. 2019, pp. 320–324.[6] 3rd Generation Partnership Project (3GPP), “Technical Specification Group Services and System Aspects; Study onpositioning use cases; Stage 1 (Release 16),” TR22.872 V16.1.0, Aug. 2019.[7] H. Wymeersch, G. Seco-Granados, G. Destino, D. Dardari, and F. Tufvesson, “5G mmWave positioning for vehicularnetworks,”
IEEE Wireless Commun. , vol. 24, no. 6, pp. 80–86, Dec. 2017.[8] A. Shahmansoori, G. E. Garcia, G. Destino, G. Seco-Granados, and H. Wymeersch, “5G position and orientation estimationthrough millimeter wave MIMO,” in
Proc. IEEE GlOEBCOM Workshops (GC Wkshps),
San Diego, CA, Dec. 2015, pp.1–6.[9] Z. Abu-Shaban, X. Zhou, T. D. Abhayapala, G. Seco-Granados, and H. Wymeersch, “Error bounds for uplink and downlink3D localization in 5G mmWave systems,”
IEEE Trans. Wireless Commun. , vol. 17, no. 8, pp. 4939–4954, Aug. 2018.[10] A. Guerra, F. Guidi, and D. Dardari, “Single-anchor localization and orientation performance limits using massive arrays:MIMO vs. beamforming,”
IEEE Trans. Wireless Commun. , vol. 17, no. 8, pp. 5241–5255, Aug. 2018.[11] R. Mendrzik, H. Wymeersch, G. Bauch, and Z. Abu-Shaban, “Harnessing NLOS components for position and orientationestimation in 5G millimeter wave MIMO,”
IEEE Trans. Wireless Commun. , vol. 18, no. 1, pp. 93–107, Jan. 2019.[12] A. Kakkavas, M. H. Castañeda García, R. A. Stirling-Gallacher, and J. A. Nossek, “Performance limits of single-anchormillimeter-wave positioning,”
IEEE Trans. Wireless Commun. , vol. 18, no. 11, pp. 5196–5210, Nov. 2019.[13] A. Shahmansoori, G. E. Garcia, G. Destino, G. Seco-Granados, and H. Wymeersch, “Position and orientation estimationthrough millimeter-wave MIMO in 5G systems,”
IEEE Trans. Wireless Commun. , vol. 17, no. 3, pp. 1822–1835, Mar.2018.[14] J. Talvitie, M. Koivisto, T. Levanen, M. Valkama, G. Destino, and H. Wymeersch, “High-accuracy joint position andorientation estimation in sparse 5G mmWave channel,” in
Proc. IEEE Int. Conf. Commun. (ICC),
Shanghai, China, May2019, pp. 1–7. [15] J. Palacios, G. Bielsa, P. Casari, and J. Widmer, “Single- and multiple-access point indoor localization for millimeter-wavenetworks,” IEEE Trans. Wireless Commun. , vol. 18, no. 3, pp. 1927–1942, Mar. 2019.[16] Y. Liu, Y. Shen, and M. Z. Win, “Single-anchor localization and synchronization of full-duplex agents,”
IEEE Trans.Commun. , vol. 67, no. 3, pp. 2355–2367, Mar. 2019.[17] T. Wang, H. Zhao, and Y. Shen, “An efficient single-anchor localization method using ultra-wide bandwidth systems,”
Applied Sciences , vol. 10, no. 1, Dec. 2019.[18] M. Rath, J. Kulmer, E. Leitinger, and K. Witrisal, “Single-anchor positioning: Multipath processing with non-coherentdirectional measurements,”
IEEE Access , vol. 8, pp. 88 115–88 132, May 2020.[19] A. Fascista, A. Coluccia, H. Wymeersch, and G. Seco-Granados, “Millimeter-wave downlink positioning with a single-antenna receiver,”
IEEE Trans. Wireless Commun. , vol. 18, no. 9, pp. 4479–4490, Sep. 2019.[20] ——, “Low-complexity accurate mmWave positioning for single-antenna users based on angle-of-departure and adaptivebeamforming,” in
Proc. IEEE Int. Conf. Acoustics, Speech and Signal Process. (ICASSP),
Barcelon, Spain, May 2020, pp.4866–4870.[21] C. Gentner, T. Jost, W. Wang, S. Zhang, A. Dammann, and U. Fiebig, “Multipath assisted positioning with simultaneouslocalization and mapping,”
IEEE Trans. Wireless Commun. , vol. 15, no. 9, pp. 6104–6117, Sep. 2016.[22] E. Leitinger, F. Meyer, F. Hlawatsch, K. Witrisal, F. Tufvesson, and M. Z. Win, “A belief propagation algorithm formultipath-based SLAM,”
IEEE Trans. Wireless Commun. , vol. 18, no. 12, pp. 5613–5629, Dec. 2019.[23] H. Wymeersch, N. Garcia, H. Kim, G. Seco-Granados, S. Kim, F. Wen, and M. Fröhle, “5G mmWave downlink vehicularpositioning,” in
Proc. IEEE Global Commun. Conf. (GLOBECOM),
Abu Dhabi, UAE, Dec. 2018, pp. 206–212.[24] J. Talvitie, T. Levanen, M. Koivisto, K. Pajukoski, M. Renfors, and M. Valkama, “Positioning of high-speed trains using5G new radio synchronization signals,” in
Proc. IEEE Wireless Commun. and Netw. Conf. (WCNC) , Apr. 2018, pp. 1–6.[25] X. Li, E. Leitinger, M. Oskarsson, K. Åström, and F. Tufvesson, “Massive MIMO-based localization and mapping exploitingphase information of multipath components,”
IEEE Trans. Wireless Commun. , vol. 18, no. 9, pp. 4254–4267, Sep. 2019.[26] R. Mendrzik, F. Meyer, G. Bauch, and M. Z. Win, “Enabling situational awareness in millimeter wave massive MIMOsystems,”
IEEE J. Sel. Areas Commun. , vol. 13, no. 5, pp. 1196–1211, Sep. 2019.[27] H. Kim, K. Granström, S. Kim, and H. Wymeersch, “Low-complexity 5G SLAM with CKF-PHD filter,” in
Proc. IEEEInt. Conf. Acoustics, Speech and Signal Process. (ICASSP),
Barcelona, Spain, May 2020, pp. 5220–5224.[28] B. Zhou, A. Liu, and V. Lau, “Successive localization and beamforming in 5G mmWave MIMO communication systems,”
IEEE Trans. Signal Process. , vol. 67, no. 6, pp. 1620–1635, Mar. 2019.[29] N. Garcia, H. Wymeersch, and D. T. M. Slock, “Optimal precoders for tracking the AoD and AoA of a mmWave path,”
IEEE Trans. Signal Process. , vol. 66, no. 21, pp. 5718–5729, Nov. 2018.[30] R. Koirala, B. Denis, D. Dardari, and B. Uguen, “Localization bound based beamforming optimization for multicarriermmWave MIMO,” in
Bremen, Germany, Oct.2017, pp. 1–6.[31] R. Koirala, B. Denis, B. Uguen, D. Dardari, and H. Wymeersch, “Localization optimal multi-user beamforming with multi-carrier mmwave mimo,” in
IEEE 29th Annual Int. Symposium Personal, Indoor and Mobile Radio Commun. (PIMRC),
Bologna, Italy, Sep. 2018, pp. 1–7.[32] H. Zhao, L. Zhang, and Y. Shen, “On the optimal beamspace design for direct localization systems,” in
Proc. IEEE Int.Conf. Commun. (ICC),
Kansas City, MO, USA, May 2018, pp. 1–6.[33] H. Zhao, N. Zhang, and Y. Shen, “Robust beamspace design for direct localization,” in
Proc. IEEE Int. Conf. Acoustics,Speech and Signal Process. (ICASSP),
Brighton, UK, May 2019, pp. 4360–4364. [34] W. W. Li, Y. Shen, Y. J. Zhang, and M. Z. Win, “Robust power allocation for energy-efficient location-aware networks,” IEEE/ACM Trans. Netw. , vol. 21, no. 6, pp. 1918–1930, Dec. 2013.[35] A. Shahmansoori, G. Seco-Granados, and H. Wymeersch, “Power allocation for OFDM wireless network localization underexpectation and robustness constraints,”
IEEE Trans. Wireless Commun. , vol. 16, no. 3, pp. 2027–2038, Mar. 2017.[36] A. Kakkavas, G. Seco-Granados, H. Wymeersch, M. H. C. Garcia, R. A. Stirling-Gallacher, and J. A. Nossek, “5G downlinkmulti-beam signal design for LOS positioning,” in
Proc. IEEE Global Commun. Conf. (GLOBECOM),
Waikoloa, HI, USA,Dec. 2019, pp. 1–6.[37] R. Vaughan and J. Bach-Anderson,
Channels, propagation and antennas for mobile communications , ser. ElectromagneticWaves. Stevenage: The Institution of Engineering and Technology, 2003.[38] T. S. Rappaport, E. Ben-Dor, J. N. Murdock, and Y. Qiao, “38 Ghz and 60 Ghz angle-dependent propagation for cellular& peer-to-peer wireless communications,” in
Proc. IEEE Int. Conf. Commun. (ICC),
Ottawa, Canada, Jun. 2012, pp.4568–4573.[39] M. T. Martinez-Ingles, D. P. Gaillot, J. Pascual-Garcia, J. M. Molina-Garcia-Pardo, M. Lienard, and J. V. Rodríguez,“Deterministic and experimental indoor mmW channel modeling,”
IEEE Antennas Wireless Propag. Lett. , vol. 13, pp.1047–1050, 2014.[40] M. Peter et al. , “Measurement campaigns and initial channel models for preferred suitable frequency ranges; DeliverableD2.1,” Mar. 2016. [Online]. Available: https://5g-mmmagic.eu/results/
SIAM Journal on Optimization , vol. 27, no. 2, pp. 616–639, Apr. 2017.[42] A. Dammann, T. Jost, R. Raulefs, M. Walter, and S. Zhang, “Optimizing waveforms for positioning in 5G,” in
Proc. IEEE17th Int. Workshop on Signal Process. Advances in Wireless Commun. (SPAWC),
Edinburgh, UK, Jul. 2016, pp. 1–5.[43] M. D. Larsen, G. Seco-Granados, and A. L. Swindlehurst, “Pilot optimization for time-delay and channel estimation inOFDM systems,” in
Proc. IEEE Int. Conf. Acoustics, Speech and Signal Process. (ICASSP),
Prague, Czech Republic, May2011, pp. 3564–3567.[44] A. Shahmansoori, G. Seco-Granados, and H. Wymeersch, “Robust power allocation for OFDM wireless networklocalization,” in
Proc. IEEE Int. Conf. Commun. Workshop (ICCW),
London, UK, Jun. 2015, pp. 718–723.[45] Y. Rockah and P. Schultheiss, “Array shape calibration using sources in unknown locations – Part I: Far-field sources,”
IEEE Trans. Acoust., Speech, Signal Process. , vol. 35, no. 3, pp. 286–299, Mar. 1987.[46] H. Messer, “The hybrid Cramér-Rao lower bound - from practice to theory,” in
Proc. Fourth IEEE Workshop on SensorArray and Multichannel Process. (SAM),
Waltham, MA, USA, Jul. 2006, pp. 304–307.[47] R. Cools, “An encyclopaedia of cubature formulas,”
Journal of Complexity , vol. 19, no. 3, pp. 445 – 453, Jun. 2003,Oberwolfach Special Issue.[48] D. Crouse, “Basic tracking using nonlinear 3D monostatic and bistatic measurements,”
IEEE Aerosp. Electron. Syst. Mag. ,vol. 29, no. 8, pp. 4–53, Aug. 2014.[49] K. Levenberg, “A method for the solution of certain non-linear problems in least squares,”
Quarterly of AppliedMathematics , vol. 2, no. 2, pp. 164–168, Jul. 1944.[50] D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,”
Journal of the Society for Industrialand Applied Mathematics , vol. 11, no. 2, pp. 431–441, Jun. 1963.[51] P. Zhang, L. Gan, S. Sun, and C. Ling, “Atomic norm denoising-based channel estimation for massive multiuser MIMOsystems,” in