Single-Anchor Positioning: Multipath Processing With Non-Coherent Directional Measurements

High-accuracy indoor radio positioning can be achieved by using (ultra) wideband (UWB) radio signals. Multiple ﬁxed anchor nodes are needed to compute the position or alternatively, specular multipath components (SMCs) extracted from radio signals can be exploited. In this work, we study a multipath-based, single-anchor positioning system that acquires directional measurements non-coherently. These non-coherent measurements can be obtained, e.g., from a single-chain mm-wave transceiver with analog beam steering or from a low-complexity ultra-wideband transceiver with switched directional antennas. The directional antennas support the separation of SMCs and the suppression of the undesired diffuse multipath component (DMC) with the beneﬁt that the required signal bandwidth can be drastically reduced. The paper analyzes the Cramér-Rao lower bound (CRLB) on the position estimation error to gain insight in the inﬂuence of the system design parameters as well as the impact of the DMC on the position error. The CRLB is compared between the non-coherent antenna setup, a conventional array with coherent processing, and a single-antenna setup. A maximum-likelihood position estimation algorithm is formulated. Its performance is evaluated with synthetically generated data as well as with UWB measurements. We show that the accuracy and robustness are signiﬁcantly improved due to the processing of angular information. Analyzing the measured data for a line-of-sight link, the median error decreases from 22 down to 7 cm, the measurements better than 20 cm increase from 46 to 95 %, and outliers above 50 cm reduce from 12 to 0 %.


I. INTRODUCTION A. MOTIVATION AND STATE OF THE ART
Fifth generation (5G) radio networks are expected to employ massive multiple-input multiple-output (MIMO) array antenna systems [1], [2] due to the dramatic increase in energy efficiency and capacity, achieved by aggressive spatial multiplexing [3], [4]. Array antennas in MIMO systems also open the gate to use millimeter-wave (mm-wave) based systems, since the array gain with a high number of antennas compensates the high pathloss present in these frequency bands  [5]- [7].
For positioning, large bandwidth enables good time resolution, whereas MIMO processing enables good angle The associate editor coordinating the review of this manuscript and approving it for publication was Vittorio Degli-Esposti . resolution by exploiting array processing [8]- [14]. This conclusion has been drawn from the analysis of the Cramér-Rao lower bound (CRLB) on the position (and orientation) estimation error(s). The CRLB is a powerful tool to analyze the performance of positioning systems with regards to the impact of system parameters, such as bandwidth, carrier frequency, and antenna configuration. Additionally, the analysis of the CRLB gives insights about the impact of model-dependent multipath component (MPC) parameters, i.e., specular multipath component (SMC) and dense multipath component (DMC) parameters [10]- [15]. An increased signal bandwidth improves the delay resolution of SMCs and also the capability to suppress the DMC [10], [12], [14], [15], and therefore it improves the delay information that can be extracted from received radio signals. Similarly, the angular information is enhanced by increasing the array aperture, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ as analyzed in detail in [11]. Here, the number of antenna elements is also related to a diversity gain in (dense) multipath channels [12]- [14], which reduces the bandwidth needed for a targeted performance. Recent works [8], [15]- [18] build on leveraging SMCs for positioning if a sufficient signal bandwidth and/or number of antennas is accessible to enable high time and angle resolution. Each SMC carries position-related information stemming from its delay only [15] or from its delay, angleof-departure (AoD), and angle-of-arrival (AoA) [16], [17], leading to a reduced need for infrastructure (number of fixed anchors) and/or an increased redundancy (robustness). However, to exploit this information, it is necessary to model the relationship between parameters that characterize the SMCs and the position/orientation of a mobile agent. A geometric model of the environment can serve for this purpose. In the most simple case, this model consists of a set of virtual anchors (VAs) -mirror images of the physical anchor at flat surfaces [8], [19] -as illustrated in Fig. 1. Multipath-based positioning systems face the challenge of acquiring such maps of environment features. The works in [16], [17] consider the impact of estimating these map features on the Fisher information of the position and orientation error. In [20]- [24], simultaneous localization and mapping algorithms are presented that estimate the environment map from consecutive measurements, achieving accurate and robust position tracking of a mobile agent in challenging indoor environments. In [8], [25], the use of mm-wave antenna arrays was discussed for multipath-based indoor positioning.
All these references show that coherent antenna array systems are advantageous for robust high-accuracy positioning.
However, their common drawback is the increased hardware complexity and the resulting (peak) power consumption. Each additional transceiver chain multiplies the hardware needed for the signal acquisition/generation and the corresponding frontend signal processing steps. This makes the use of parallel, coherent processing of all the antenna signals very expensive, in particular for ultra-wide bandwidth (UWB) signals with sampling frequencies in the GHz-range. Hybrid analog-digital architectures have been proposed as an alternative for massive-MIMO mm-wave base-stations aiming at very high antenna directivities [26]- [28]. In this architecture, the number of radio frequency chains is much smaller than the number of antennas [29]; potentially only one transceiver chain is implemented. The main idea is to use analog beamforming (e.g. phase shifters or directive antenna elements) to gain directivity and only process one or a few received signals coherently [30].
A similar principle has been employed in [31], where a single-anchor, multipath-based positioning system has been described for high-accuracy (indoor) positioning for low-cost Internet-of-Things applications. A switched antenna system with four directive antenna elements is used, pointing in the four cardinal directions, in order to facilitate the reliable separation of MPCs [32]. As an example, the top left of Fig. 1 illustrates the antenna configuration from [31], [32], where UWB radios are used according to the IEEE802.15.4 standard.
With the newly proposed ''phase-based ranging'' extension [33], accurate delay and angle estimation also becomes applicable to Bluetooth Low Energy signals. However, it remains to be investigated if a multipath-based positioning system can be realized with this technology. Furthermore, switched directive antennas were also proposed for RSS-based indoor positioning [34], [35]. In [36], the CRLB was derived for RSS-based indoor positioning using switched directive antennas.
Both analog/hybrid beamforming as well as switched antenna systems are not capable of acquiring coherent measurements for conventional phase-based array processing, because signals from multiple antennas are measured consecutively rather than in parallel. Phase coherence between such consecutive measurements is very hard to achieve, given the high demands on radio-frequency oscillators and the potential dynamics of the environment, especially with mm-wave radios [28].

B. CONTRIBUTIONS OF THE PAPER
In this work, we analyze a single-anchor positioning system exploiting SMCs, which is capable of acquiring channel impulse response measurements with a set of different directive antennas. Potentially, these measurements are not phase-coherent with respect to one another. We derive the CRLB on the position error to evaluate the theoretical performance limit of this positioning system, in comparison to a conventional antenna-array system and a single-antenna setup. An illustration of both antenna setups is shown in Fig. 1. For conventional coherent processing, the anchor is equipped with an array of antennas of known geometry radiating in an isotropic manner (shown on the right for a rectangular constellation). For so-called ''non-coherent'' processing, a set of directional antennas is used with known beampatterns (as shown on the left). We also consider the DMC to model the interference effect of non-specular multipath components in a more realistic form in contrast to an AWGN model.
We significantly extend our initial error bound analysis provided in [32], in particular, a detailed analysis of AoA information has been included yielding guidelines for the design of antenna radiation patterns. Furthermore, we describe a positioning algorithm based on our previous work in [31], [37], [38] and analyze its performance in comparison to the CRLB, using synthetically generated data as well as measured data. Our specific contributions are: • We analyze the CRLB for a multipath signal model, considering non-coherent antenna arrays with directional beampatterns in comparison to conventional coherently processed antenna arrays (Sections III & IV).
• We quantify the contributions of delay and angle information of SMCs to the position information, taking into account self-interference by the DMC (Section V).
• We develop and analyze positioning algorithms for the non-coherent directional antenna array (Section VI).

C. NOTATIONS
Boldface upper case letters represent matrices. Boldface lower case letters denote column vectors. Superscripts T , * and H denote matrix transpose, complex conjugation and Hermitian transpose, respectively. The Kronecker product is denoted with ⊗. · is the Euclidean norm. | · | represents the absolute value.Â denotes an estimate of A. I [ · ] represents the identity matrix with dimension denoted in the subscript [ · ]. E · denotes the expectation operator.

II. PROBLEM FORMULATION
We consider the task of finding the position p of an agent node using radio signal measurements from one anchor node located at known position a 1 . We first examine the complete channel model and extract an approximation that contains position-related parameters. Then, we consider the received signal model to describe the recorded observations that will be used to determine the agent position. Finally, we describe how the channel parameters are related to the geometry of the environment.

A. CHANNEL MODEL
The radio channel between the agent at position p and the anchor at position a 1 is described with the spreading function [39] according to where δ(·) denotes the Dirac delta function. This equation describes the superposition of MPCs that originate from reflections in the environment. For the sake of simplicity, we assume a two-dimensional scenario with horizontal-only propagation (in the azimuth plane). 1 The rational behind the selection of the model in (1) is that for radio measurements conducted with a finite observation aperture in space and frequency, not all MPCs can be resolved in distinct SMCs. Therefore, we assume that separated MPCs are collected into a set of k = 1 . . . K SMCs and all unresolvable MPCs originating for example from diffuse scattering are described by a DMC [39], [40]. Each SMC is described by amplitude α k , phase ζ k , AoA φ k and delay τ k . The latter two parameters can be related to the agent position via the geometry of the environment as exemplified in Fig. 1 (for the delay, the corresponding path length d k is shown). This will be described in detail in Sec. II-D. The DMC ν(φ, τ ) ∈ C is modeled as a complex circular (i.e. zero-mean) Gaussian random process [40], [41]. Assuming uncorrelated scattering (US) in the delay and angular domains, the auto-correlation function of ν(φ, τ ) is given by where S ν (φ, τ ) describes the azimuth-delay power spectrum [39] at the anchor position.

B. RECEIVED SIGNAL MODEL
The anchor node employs an antenna array which consists of M antennas, each of which exhibits a beampattern b m (f , φ). When referring to the general anchor position, we use a 1 which is defined to be the mass point of the antenna array.
On the other hand, the position of the mth antenna element at the anchor is denoted by a (m) 1 . The agent transmits a lowpass-equivalent signal s(t) modulated by carrier frequency f c and the anchor receives signal r m (t) using antenna m. 2 We aim for a compact description and introduce the following assumptions.
• For the antennas at the anchor, in frequency domain, we assume identical beampatterns over all relevant frequencies, i.e., we use b m (f , φ) = b m (φ) for the complex-valued azimuth antenna gains.
• For the time-domain description of radio waves impinging at the antenna elements, we make use of the far-field plane-wave assumption, i.e., we assume planar wave fronts instead of spherical ones.
• Introduced time delays at antenna elements in relation to the mass point will be only considered in terms of a phase change, i.e., the time delay of the signal envelope will be neglected. This is usually referred to as 1 An extension to three dimensional scenarios with horizontal and vertical propagation is straightforward, but it would lead to cumbersome notations without bringing significant additional insights. 2 Note that we assume that the anchor and the agent are synchronized. This can be achieved by using a two-way transmission protocol [42]. However, the proposed model can be extended to non-synchronized anchor-agent links along the lines of [15] (based on the fact that the relevant geometric information is also contained in the time differences of the SMCs).
the narrowband/wideband assumption 3 and it is fair to use since the envelope information is negligible with respect to the phase information, especially in practical situations [11], [13], [14]. We apply these assumptions to obtain the received signal at antenna m via convolution of the spreading function in (1) with the transmitted signal s(t) and antenna response b m (φ)e jζ (m) (φ) resulting in which can be separated in three distinct parts: The first part contains the position-related SMC parameters {φ k } and {τ k } which shift and scale the transmitted signal. The second part r DM m (t) is a stochastic process characterizing the self-interference due to the DMC. Finally, measurement noise w m (t) is modeled as additive white Gaussian noise (AWGN) with double-sided power spectral density of N 0 2 . The SMC-related phase shifts are given as where ζ (m) (φ) captures the antenna related phase shift and ζ k the radio channel related one.
In the case of coherent array processing, only the K radio channel related phase shifts ζ k are unknown in (4) and the systematic phase-offset introduced by the placement of antenna m w.r.t. reference position a 1 is given as where τ (m) (φ) is the time delay due to the distance of the antenna elements (situated at a (m) 1 ) relative to the phase center at a 1 , d (m) and φ (m) denote the distance and angle of the mth antenna w.r.t. a 1 (cf. Fig. 1), and λ is the wavelength at f c . In the case of non-coherent processing, all measurements at the M antennas have unknown phase-offsets meaning that all K · M phase shifts ζ k,m given by (4) are unknown.
The second part of (3) characterizes self-interference caused by the DMC and it is given as We look at correlation properties of this signal for a pair of antennas where we applied the US assumption from (2). It should be noted that this assumption also allows us to use the same stochastic process ν(φ, τ ) for all antenna elements, because the assumption implies homogeneity in the spatial domain [44,Ch. 2.4].
Here, we have on the one hand the SMC-related term x(θ) as a function of the SMC parameters θ , and on the other hand the DMC and noise related term n. In the following, we will describe these two terms in detail.

1) SMCs
The SMC part of the observation in (7) is described by The parameter vector θ contains the SMC parameters where its components are of dimension R K ×1 . Of special mention is the phase parameter ζ , which is only of the same dimension if coherent processing is performed, i.e., if the antenna phases are determined by (5). However, our focus will be on the non-coherent case where ζ ∈ R (KM )×1 contains all phase shifts ζ k,m , since they are all unknowns. Hence, in this case, the stacked parameter vector θ is of dimension R (3K +KM )×1 .

2) DMC AND NOISE
The vector n = n ν + w ∈ C MN ×1 represents the DMC process and the measurement noise as a Gaussian process with covariance matrices where σ 2 w = N 0 /T s is the measurement noise variance and C v is the DMC covariance matrix. For the DMC covariance matrix we need to take correlations between antenna elements into account. To this end, we can separate C v into blocks of size N × N , each describing the correlations between one pair of antennas indexed as, e.g., (m, m ). We may use (6) to write one block of the DMC covariance matrix as

D. RELATING SMCs TO GEOMETRY
We describe the SMC parameters contained in θ in more detail, considering their relation to the agent position. It should be emphasized that we will use the environment geometry, that is, we assume that each SMC involves the reflection on flat surfaces such as wall segments, as shown in Fig. 1. These segments can be described by a surface normal, or rather in the case of azimuth plane operation, by a single segment angle as will be described later. We aim for a compact and efficient computation of the parameters when anchor and surface positions are known. The SMC delay (time-of-flight) τ k is related to the path length d k by where c is the propagation velocity. Here, a k is the position of a VA [8], [45] that results from mirroring the anchor position on a known reflective surface in the environment, as exemplified in Fig. 1 by a 2 and a 3 . Note that Fig. 7.10 in [19,Chapter 7.5] provides a detailed description of this operation called the image-source principle. For higher-order reflections, the mirroring process is repeated in the sequence of reflecting surfaces along the propagation path.
Considering the angle domain, the AoD ϕ k at the agent is given by ϕ k = (a k − p). Using the angles of all involved wall segments, the AoD ϕ k at the agent can be related to the AoA at the anchor node by where O(k) denotes the order of the reflection associated with MPC k, and φ (k) seg j denotes the angle of the jth involved reflective segment in the propagation path with j = 1, . . . , O(k) ordered according to the sequence of bounced surfaces. As an example, see Fig. 1 for the segment angle related to VA a 3 .
What remains is the amplitude α k and phase ζ k,m . The amplitude incorporates effects such as path-loss and reflection losses, whereas the phase incorporates the influence of reflections depending on the involved materials. These parameters are not easily related to the geometry and are thus treated as nuisance parameters in the position estimation problem.

III. FISHER INFORMATION OF CHANNEL PARAMETERS
On the basis of the signal model in (7) and (8), we can examine the useful information present in observed signals. To this end, we will determine the Fisher information matrix (FIM). The elements of this matrix quantify the amount of information that the observable vector r carries about the unknown parameters θ .

A. FISHER INFORMATION MATRIX
The FIM I θ on the parameter vector θ is used to obtain the Cramér-Rao lower bound (CRLB) [46], [47] for an estimated parameter vector ofθ via The CRLB provides a lower bound on the variance of any unbiased estimator for the model parameters. The FIM, considering the Gaussian model with known covariance matrix C n , is defined as [47,Chapter 15.7] where J θ denotes the Jacobian of the signal model with respect to the elements of a length-L θ parameter vector In our case, using non-coherent processing, the number of parameters is L θ = K (3 + M ). We define sub-matrices I ϑϑ [15], where ϑ, ϑ ∈ {φ, τ , α, ζ }. These sub-matrices cover all combinations of parameter types to assemble the full FIM from (15). The structure of the full FIM I θ using these sub-matrices can be found in (49), whereas one element [·] k,k of I ϑϑ is outlined in (50), both shown in Appendix A.

B. TOWARDS CRLB: EQUIVALENT FIM
Our main interest concerns the lower bound on the estimation of the delays τ k and AoAs φ k , and the resulting position estimation performance. Hence, we want to determine sub-matrices of the inverse FIM, in particular the respective block matrices on the main diagonal related to one type of parameter. For this purpose, it is beneficial to define the equivalent FIM (EFIM) [10] as where ϑ stacks the vectors of all remaining parameters, i.e. ϑ = θ \ ϑ. The CRLB on the corresponding parameter can be obtained by computing the inverse of the respective EFIM. It should be noted that we rearrange the sub-matrices in (49) such that any desired parameter type is in the top left to make use of this inversion lemma. Numeric evaluation of (17) will allow for the evaluation of the CRLB for the respective parameter vector. The benefit of the EFIM definition is that the cross-dependence of the estimation of the parameters in ϑ on any other parameters, described by the matrix I ϑϑ , turns approximately to 0 in many cases. Thus, it becomes sufficient to derive the inverse of the sub-matrix I ϑϑ , yielding far more insightful formulations. In the following, we will present the derivations of the EFIM for the SMC delay and angle estimation using the assumptions given in Appendix B that allow further insightful discussions in terms of contributions to position information. Note that in general (17) can be used to compute the CRLB on the position error numerically without using the assumptions given in (64), (65) and (68). However, to gain more insights into the structure of the CRLB, we use these assumptions in what follows.

C. EFIM FOR SMC DELAY ESTIMATION
Based on the assumptions (64), (65) and (68) in Appendix B, the EFIM for delay estimation can be derived as (61) given in Appendix A. It can be rearranged in the following form is the (mean-square) bandwidth of the whitened signal s(τ k ) andṡ(τ ) = ∂s/∂τ denotes the derivative of the signal s(τ ) w.r.t. the delay τ , and is the signal-to-interference-plus-noise-ratio (SINR), quantifying the level of DMC plus noise interference in one SMC. We use · 2 H to denote the squared weighted norm (taking the noise model into account) which is described by (67) in Appendix B. The factor ξ delay k ∈ [0, 1] can be interpreted as the information loss on the SMC delay estimation induced by the complex amplitude estimation in the DMC [12], [15], and it is given in (61) in Appendix A. We introduce the whiten- is the mean-square bandwidth of signal s(τ ), and the effective SINR, SINR k = SINR k γ k ξ delay k . Then, the EFIM for delay estimation in (18) can be rewritten as which is equivalent to what was shown in [12], with the difference that there is a weighting by the squared beampattern values, hence a gain dependent on the directionality.

D. EFIM FOR SMC AoA ESTIMATION
We follow an analogous approach for the EFIM of the AoA estimation, i.e., we use the assumptions (64), (65) and (68) in Appendix B. The EFIM for AoA estimation in (62), in Appendix A, can be rearranged in the following form can be interpreted as the information loss on the SMC AoA estimation induced by the complex amplitude estimation in the DMC, and it is given in (62) in Appendix A.
A common property of antenna characteristics is an approximately symmetric beampattern, i.e., the beampattern can be described by an even function. In this case, due to the property of even functions exhibiting an odd function as derivative, there is orthogonality between the function and its derivative amounting to ḃ (φ k ), b(φ k ) ≈ 0. Therefore, we have ξ angle k ≈ 1 in (62) neglecting the information loss related to the angles in the upcoming derivations. Introducing the non-coherent normalized square array aperture the EFIM for AoA estimation in (22) can be rewritten as

E. COMPARISON TO COHERENT PROCESSING
For a conventional antenna array with coherent processing, we obtain (cf. [11], [13]) with which will serve as a reference result. The factor (26) has been interpreted as a coherent normalized squared aperture that scales the amount of angle information available. (Note the similarity to the scaling of delay information by the squared bandwidth.) The accompanying array response beampattern using (5) is given by For a numeric comparison of the squared apertures from (23) and (26), it is convenient to write the beampattern of antenna m as a Fourier series, where coefficients c describe the generic beam pattern b(φ) which is rotated by m 2π M for antenna m. For non-coherent processing, a real-valued beampattern is assumed, hence where the inner sum is the autocorrelation of the Fourier coefficients, c , at lags ηM . For η = 0, this expression corresponds to a constant squared norm, whereas for |η| ≥ 1, deviations from this constant are quantified. The deviations have a periodicity of 2π M , obviously. We can use this insight to design antenna patterns that minimize or avoid fluctuations of b(φ) 2 as a function of φ, which is important to yield uniform delay information, independent of φ k , cf. (18). Such constant delay information will be obtained when only coefficients c = 0 for | | ≤ (M − 1)/2, yielding b(φ) 2 = M |c | 2 . Normalizing the energy of the beam pattern to one, i.e. b(φ) 2 = |c | 2 = 1, we get b(φ) 2 = M . It is seen that an SINR gain by a factor of M is obtained, which is identical to the SINR gain obtained with conventional coherent processing [12], [13].
To gain insight in the equivalent aperture value of the non-coherent antenna configuration, we analyze the expression (23). We re-use the Fourier series representation of b m (φ) from (28), to reformulate (23) as where the inner sum is the autocorrelation of the Fourier coefficients of the beampattern derivative, c , at lags ηM . Again, this expression corresponds to a constant for η = 0, whereas for |η| ≥ 1, deviations from this constant are quantified, which are periodic with 2π M . We can use this insight to design antenna patterns that minimize or avoid fluctuations of the angle information as a function of φ. In particular, constant angle information will be obtained when only coefficients c = 0 for | | ≤ (M − 1)/2, yielding Choice of a constant c = 1/ √ M for | | ≤ (M − 1)/2 yields a Dirichlet kernel for the beam pattern with normalized energy b(φ) = 1. It will have an equivalent squared aperture which scales quadratically in the number of array elements.  Fig. 2b we show a choice of coefficients that results in a raised-cosine shape, which means reduced sidelobes for the beampattern shape at the cost of a ripple in the equivalent aperture, i.e., varying delay and AoA information for different directions.
Comparison to a uniform linear array with coherent processing and λ/2-spacing [13, Eq. (33)] reveals the same scaling with M and an increase by a factor of π 2 in broadside direction. The improvement is apparently at the cost of information in end-fire direction. In a uniform circular array, a constant D 2 λ (φ) ≈ M 2 /(16π 2 ) is obtained (for M 1) which differs from (32) by a factor of 3. Fig. 3 shows illustrations of these array geometries with accompanying array responses (see [48,Sec.3]) and information gains for a steering in direction π 2 . These results demonstrate the close correspondence between the angle information retrieved from the non-coherent antenna setup and a conventional, coherent array. The root-EFIM of the angle parameter will be reduced  by a factor of √ 3 for the non-coherent case. We will further evaluate and compare the two antenna configurations in Section V by means of numerical results.

IV. ERROR BOUND FOR POSITIONING
This far, we derived the full FIM and EFIMs for our parameter vector θ . In the following, we evaluate the obtainable position accuracy for our measurement model. For that VOLUME 8, 2020 matter, we determine the CRLB for the agent position, termed position error bound (PEB) [10]. The PEB acts as a comprehensible quality measure that can be used as a benchmark for practical estimators.

A. GENERAL PEB
When the FIM is available for the observation model parameters, we can relate it to the PEB via a parameter transformation [10]. As a first step, we introduce a new parameter vector Here, we have the actual desired quantity, the agent position p, but we still have to include the SMC amplitudes and phases as nuisance parameters. Since we have D = 2 dimensional positions and again phases ζ k,m , we obtain L ψ = D + K (M + 1) for the respective number of parameters. We relate the original parameter vector to the new one via the Jacobian where matrix J p describes the spatial gradients defined by Together with the parameter FIM from (49), we obtain the FIM for the new parameters The CRLB for ψ is related to the inverse of this matrix (cf. (14)). However, we are interested in the positions only, hence we look at the trace of the top left of the inverse that defines the PEB with Taking the square of the right-hand side results in what is commonly known as the squared PEB, which quantifies the lower bound on the variance of the absolute position error.

B. SPATIAL GRADIENTS
For the spatial gradients (35), we apply the derivative with respect to the agent position to the parameter definitions in (12) and (13), resulting in ∂ ∂p (α k ) = 0 ; ∂ ∂p (ζ k,m ) = 0 where e(ϕ) denotes a unit vector pointing in the direction of the AoD ϕ. Furthermore, we assume that the amplitudes α do not depend on the agent position p. Since the phases ζ stem from reflection properties and the measurement device, it is evident that they are also unrelated to the agent position. With this, we can assemble the Jacobian from (34) and follow the derivations until (37) to obtain the PEB.

C. PEB USING EFIM
We use the EFIMs (21) and (24) to get a closed-form solution for (37), in order to obtain an insightful formulation. To make use of the EFIMs, it is beneficial to define subsets of the Jacobian from (35) that contain only the position-related parameters described by The parameter transformation (36) enables the use of the inversion lemma (17) to compactly obtain the EFIM for the agent position via With this, the full EFIM for the agent position is described by where D r (ϕ) = e(ϕ)e T (ϕ) is the ranging direction matrix [10]. For each SMC, we can see an information gain in radial direction via delay information (quantified by the squared bandwidth and the speed of light, i.e., β 2 k /c 2 ) and in tangential direction via angle information (quantified by the non-coherent squared antenna aperture and the SMC distance, i.e., D 2 b (φ k )/d 2 k ).

V. NUMERICAL RESULTS
This section aims at a numeric validation of the potential of the beampattern-enhanced, non-coherent positioning system. We evaluate the PEB for specific scenarios in indoor environments. Three setups are used to compare the performance obtained by the proposed non-coherent array with a conventional coherent array and a reference setup with a single antenna.

A. EVALUATION SETUP
We assume given wall segments that form rooms with simple geometries and a given anchor position a 1 . The number of dedicated wall segments is denoted by K seg . The anchor position is used together with the wall segments to determine VAs via mirroring operations, as described in Section II-D. With the VA positions, the parameters τ k and φ k are determined as described by (12) and (13), respectively. We consider up to second-order reflections, which means the number of considered SMCs is given by K ≤ K seg (K seg − 1) + 1. Note that for each agent position, VA visibility tests are performed, which may reduce this number on an individual basis [49]. The real-valued amplitudes α k are simulated using the free space path-loss model, whereas each reflection results in a loss of 3 dB. Additionally, phases are applied to the amplitudes according to e j2πf c τ k . For the agent positions, a uniformly sampled grid is spanned over the floorplan of the room with a spacing of 2 cm. To provide a context to previous work, we consider an ''L-shaped'' room, as shown in Fig. 5, which was also used in [15]. Hence, with K seg = 6 wall segments present, the number of usable SMCs amounts to K ≤ 31. We also consider the same pulse signal for s(t) as in [15], using a rootraised-cosine waveform with roll-off β = 0.6 and bandwidth 1/T p = 1 GHz at a carrier frequency of f c = 7 GHz. For the DMC, we use the assumptions in (64) and (68), resulting in a block diagonal covariance matrix C n with identical blocks for each antenna. We use a delay power spectrum S ν (τ ) which exhibits a double-exponential shape as defined in [41] by Eq. (9) and illustrated in Figs. 4a and 4b. The used parameters are γ 1 = 20 ns, γ rise = 5 ns and χ = 0.98, representing the decay and rise exponents and the relative power at excess delay zero, respectively. The power parameter 1 is chosen such that we reach a Ricean K-factor of |α To define the signal-to-noise ratio (SNR) of the synthetically generated measurements, we introduce the SNR at 1 m, i.e., which relates the received LOS component energy at a distance of 1 m from the anchor to the noise power spectral density. Additionally, we define the receiver SNR as where the parameters α k and τ k are determined for a specific agent (''receiver'') position. This quantity takes distance and visibility conditions of the SMCs into account (e.g., K might be reduced due to non-visible SMCs) and thus relates to the impact of the propagation channel for a selected measurement noise level. Fig. 4 illustrates the channel parameters (SMC and DMC) for two example agent positions, in Fig. 4a for an agent position where the LOS and four SMCs are visible, and in Fig. 4b for an agent position where the LOS is blocked and only two SMCs are visible. For visualization purposes, only first-order SMCs are considered, i.e., we have K ≤ 7. The AWGN level is shown for two N 0 values where SNR (1m) = 30 dB and SNR (1m) = 50 dB. In Sec. VI, a positioning algorithm will be presented and evaluated on these simulated data.

B. PEB RESULTS FOR DIFFERENT SETUPS
We evaluate the PEB described in (37) for each agent position, where we use both the full FIM, whose elements were computed using (49) (the respective sub-matrices are given in (51)-(56)), as well as the simplified, canonical EFIM from (38). We also select representative agent positions to show the respective scaled error ellipses that illustrate the components of the PEB in 2D-space. Throughout this section, we choose N 0 (and in turn σ 2 w ) such that we get SNR (1m) = 29.5 dB as a reference value. Setup 1: Single omni. First, we recreate the results from [15], examining the PEB for the multipath model considering only a single omni-directional antenna at the anchor. Hence, we set M = 1 and b m (φ) = 1 √ 2π . The resulting PEB values are shown in Fig. 5a. While the setup allows for overall low PEB values due to the good resolution in delay domain, we can clearly see regions where the bound sharply decreases, caused by non-resolvable SMC path overlap. Fig. 8 shows the cumulative frequency of PEB values over the agent positions to quantitatively compare the positioning performance for the results illustrated in Fig. 5. The degradation in positioning accuracy due to path overlap is visible when considering the 90 percentile and above, where the error increases to multiple decimeter.
Setup 2: Array processing (coherent). We increase the antenna number to M = 4 and assume a conventional array with coherent phase processing. The antennas are on a circle with a constant radius d (m) and φ (m) = π 2 · m, such that an inter-antenna spacing of λ/2 is achieved, as illustrated in Fig. 3b. To use phase coherent processing, this setup requires the antenna phases ζ m (φ) to be determined by the respective AoA φ. Hence, for this setup, we change the vector ζ to only include the K random SMC phases ζ k , cf. (1). VOLUME 8, 2020 The PEB for this constellation is shown in Fig. 5b. This time, the angle information helps resolving SMC path overlap while also increasing the overall accuracy, leading to a PEB that is dominated by SMC visibility properties.
Non-resolvable path overlap remains only in a few regions, due to a smaller set of visible SMCs.
These performance improvements compared to the single omni setup are quantified in the cumulative frequency plot shown in Fig. 8, specifically when regarding the 90 percentile where sub-decimeter levels are achieved.
Setup 3: Directional beampatterns (non-coherent). Finally, we apply the proposed non-coherent array of M = 4 directional antennas. We use a Dirichlet kernel as defined in (28), described with M non-zero Fourier coefficients c = 1/ √ M for ∈ {−1.5, −0.5, 0.5, 1.5}, which is illustrated in Fig. 2a. Each antenna uses the same pattern but is pointed towards one of the four cardinal directions. For this setup we consider the full vector of K · M unknown phases ζ k,m . The resulting PEB is shown in Fig. 5c. The achieved accuracy comes close to Setup 2, i.e., directional antennas show a similar capability of resolving the path overlap, only slightly worse than the coherent omni-array setup. This observation is also confirmed when regarding the quantitative evaluation shown in Fig. 8, where we note the resemblance between the results of the coherent and non-coherent setups.
Additionally, in Fig. 5d, the PEB is shown for the simplified, canonical EFIM (38). It can be seen that the simplified PEB follows closely the full form, justifying the assumption in (65). The only significant differences are present in selected regions where SMC path overlap occurs. For the top left position, it is shown how the Fisher information from (38) affects both the radial and tangential position error via the information gains from delay and angle information, respectively.
Time vs. angle information. With the good achievable accuracy of the non-coherent setup, we want to determine the main contributor to the bounds in terms of position-related parameters. Specifically, we aim to show how much information the SMC delays add compared to the SMC angles. Also, we are interested in the effect of an increased number of antennas at the anchor, which puts the system in the context of mm-wave pencil-beam setups. We use the full FIM of the non-coherent model with M = 4 and M = 16 antennas to evaluate the PEB using only delays or angles, respectively, where the respective unused parameter (either τ k or φ k ) is treated as another nuisance parameter.
The PEB results, considering delay information only, are shown in Fig. 6. From Fig. 6a, it is evident that the SMC delays carry the main contribution to the achieved accuracy, because the results come close to what was achieved with the full parameter set (cf. Fig 5c). Fig. 6b shows the same evaluation using M = 16 antennas, which achieves overall slightly better results, however the improvement is not significant.
In contrast, Fig. 7 shows the PEB results considering angle information only. Please note that the PEB is now shown in meter, due to the significanty higher error values compared to the delay information case. The error ellipses show that, predominantly, tangential information is provided. Furthermore, a large degradation is evident with an increased distance from the anchor. This is a well-known drawback of angle measurements, which can be seen mathematically from the inverse distance-scaling of the angle term in (38). Using M = 16 antennas greatly increases the positioning performance, which stands in stark contrast to the delay information improvements. The 3-fold standard deviation ellipses show that, for positions in LOS conditions, the tangential deviation is minimized. Note that, for visualization purposes, these ellipses are scaled by a larger factor for M = 16, due to the significantly lower PEB values. Fig. 9 shows again a quantitative evaluation in terms of cumulative frequency for the PEB values of the described setups using delay and angle information respectively (cf. Figs. 6 and 7). With delay information, a significantly better performance is achieved. This confirms the observations from the qualitative PEB analysis, which identifies the delay information as the main contributor to accurate position estimates, justified by the large bandwidth. Increasing the number of used antennas to M = 16 decreases the error consistently. For the delay information case the accuracy scales by a factor of about 2. Considering the angle information case, the number of antennas has a higher impact. Here, using M = 16, the error decreases by a factor of 10. It is expected that the angle information will surpass the delay information when M is increased by a factor of 8 to M = 128. However, this result also depends on the room geometry, since angle information drops with d 2 k .

VI. POSITIONING ALGORITHM
In this section, we leverage the gained insights and present a positioning algorithm capable of estimating the agent position p given a measurement r. We derive the algorithm based on the received signal model (7) where we identify τ , φ, α and ζ as unknown variables. While τ and φ can be expressed as function of p (see Section II-D), we need to estimate the nuisance parametersα ≡α(r, p) andζ ≡ζ (r, p) jointly withp. The positioning is performed by a maximum likelihood (ML) approach, on the basis of the log likelihood function L(r|p, α, ζ ) derived from (7), which will be formulated later. Two algorithm variants will be given that vary in the noise model. The ML optimization problem can be formulated asp where the estimates (α,ζ ) result from evaluated for one specific positionp. The non-linear relations in (41) prevent a closed-form solution forp. Moreover, an independent estimation ofα,ζ is infeasible which sacrifices computational efficiency. As a remedy, we propose to evaluate the log likelihood function in (41) VOLUME 8, 2020  numerically for positionsp ∈ P within the communication range. A straight-forward choice for the set of evaluation positions P is a reasonably dense position grid spanned over the floorplan of the considered environment, whereas the granularity of this grid poses a lower limit on the position error. Alternatively, prior range estimates can be used to reduce the set to a sphere around the anchor as proposed in [37] and evaluated in [31]. Then, to obtain the position estimate via (41), the likelihood is evaluated for all the chosen agent positions. To further enhance the feasibility of the algorithm, we relax the joint estimation in (42), as shown in the following.

A. DERIVATION OF THE POSITIONING ALGORITHM
An independent estimation of α and ζ is intractable for the original optimization problem, but possible by assuming that the measurements r m are independent across m and by assuming the absence of overlapping SMCs. Following the derivations from [38], we start by factorizing L(r|p, α, ζ ) using assumption (63), resulting in with L m (r m |p, α, ζ ) = −det{C n The covariance C n (m) of the mth antenna is defined in (64) and the vector x m (p, α, ζ ) ∈ C N is where we introduced an auxiliary variable α k,m = α k e jζ k,m .
For an independent estimate of the parameters in α and ζ , we assume that specular reflections do not overlap in the delay domain as stated in assumption (65). Then, maximizing (44) Subsequently, we maximize (43) w.r.t.ζ k,m , contained inζ , and the results can be expressed aŝ usingα k,m obtained from (46). The SMC amplitudesα k result from maximizing (43) w.r.t α k . Again, we assume no SMC path overlap, and obtain after some derivationŝ The steps performed by the position estimator can be summarized as follows. To evaluate L(r|p,α,ζ ) from (41) for a specificp ∈ P, first, SMC parameters τ and φ are calculated using (12) and (13) respectively, followed by estimatingα k,m via (46). Subsequently,ζ andα result from (47) and (48), respectively. Finally, these estimates are plugged into (43) to calculate L(r|p,α,ζ ). This procedure is repeated for all p ∈ P andp with the highest associated likelihood is chosen as position estimatep (cf. (41)).

B. EVALUATION OF THE POSITIONING ALGORITHM
We assess the performance of the presented positioning algorithm by comparing the mean-squared error (MSE) between the true p and estimated positionp i , on the basis of observation data taken from N MC realizations of r m , where we use both synthetically created data and measurement data. We compare the performance using two kinds of estimators which vary depending on the used noised model: First, we follow the derivations in Section VI-A where we assumed knowledge regarding the noise covariance, described by C n (m) (see (44)), resulting in an algorithm that makes use of the DMC statistics, hence we refer to this algorithm as DMC-based (note that the AWGN is also included in the covariance matrix). 4 Second, as a simpler alternative, we use an estimator with non-accessible C n (m) where no DMC is used (this follows the algorithm described in [38]), which can be accomplished by replacing C n (m) in (44) and (46) with an identity matrix I N . We refer to this variant of the algorithm as AWGN-based.

1) SYNTHETICALLY GENERATED DATA
In the first evaluation, we obtain observations by performing N MC = 1000 Monte-Carlo runs. We focus on the potential to estimate positions in LOS as well as in non-LOS conditions, located in the L-shaped room from Section V where agent positions [6,6] T and [3.7, 1.8] T correspond to the LOS and non-LOS case. The signal parameters and beampatterns are identical to Section V. For the SMCs, we use only first-order reflections, hence, we have K ≤ 7. The observations are created using (7). In the first experiment, we keep the DMC power constant. The experiment is performed for varying AWGN levels, whereas the noise variance σ 2 w = N 0 /T s is set using (39) such that SNR (1m) values in the range of 10 to 60 dB are obtained. Fig. 4a illustrates the used SMCs in comparison with DMC and AWGN for the two agent positions and two AWGN levels. Fig. 10 presents the MSE achieved comparing the DMC-based (solid blue) and AWGN-based (dashed cyan) estimators using M = 3 antennas for various AWGN levels. The corresponding PEB (black) from (37) is shown for comparison. In general, we can observe that the PEB for the LOS case is lower, showing the vital position information contained in the LOS component. In non-LOS conditions the achievable accuracy is decreased and moreover, the PEB is approached at a higher SNR (1m) value of 45 dB (SNR r = 24 dB) in comparison to LOS conditions where the PEB is approached at SNR (1m) = 25 dB (SNR r = 15 dB). A comparison between the two proposed estimators illustrates the importance of considering the DMC power in the estimation procedure. At SNR (1m) values above 50 dB, the position error of the AWGN-based estimator saturates  approach similar values. Increasing SNR (1m) (and thus SNR r ), both PEBs decrease but PEB p saturates due to the constant DMC power. This saturation demonstrates again that dense multipath effectively limits the achievable accuracy regardless of the AWGN levels of the observations. It has to be emphasized that this unavoidable effect of the radio channel itself cannot be reduced by increasing SNR (1m) . However, the insight from the theoretical results shows how it can be eased by increasing the signal bandwidth and the number of the antennas or by using directional antennas.
In the second experiment, we evaluate the impact of number of antennas M . The beampattern is obtained using (28) with M non-zero Fourier-coefficients c = 1/ Fig. 11 illustrates the MSE for LOS and non-LOS conditions achieved by the DMC-based (solid blue) and AWGN-based (dashed cyan) estimators. For the AWGN, we set N 0 such that SNR (1m) = 40 dB, whereas the remaining parameters are equal to the first experiment. In general, we can observe that the PEB is less sensitive to M , e.g. increasing M by a factor of four accompanies with a reduction of the PEB by a factor of two, as anticipated in Section III. However, at large M the position estimator has more independent observations which enhances the suppression of the DMC, as demonstrated in the non-LOS case.

2) MEASURED DATA
The second evaluation is performed on the basis of a measurement campaign described in [32] and [38], which was conducted in a 6 × 8 m laboratory room with a single anchor at VOLUME 8, 2020 position a 1 and an agent at position p as illustrated in Fig. 12. We consider again first-order reflections only, hence, for this room we have K ≤ 5. Four directive antennas were used at the anchor, whereas the complex-valued beampatterns b m (φ) were available as a codebook with a resolution of 10 • . Measurements to the agent were performed using an Ilmsens Correlative Channel Sounder [51]. Each measurement was convolved by a raised cosine pulse with a pulse width of T p = 2.4 ns, a roll-off factor of β = 0.9, and a carrier frequency of f c = 5.4 GHz. The agent was placed at positions p on a 15 × 14 grid with 5 cm spacing as shown in Fig. 12, resulting in 210 measurements in (unobstructed) LOS condition. Both, anchor and agent, were placed at a height of 1.5 m, leading to reduced floor and ceiling reflections, since all involved antennas exhibit narrow elevation patterns. For the DMC-based algorithm, we estimatedĈ n (m) individually for each antenna using the 210 measurements to estimate the DPSŜ ν (m) (τ ). For each grid position, the algorithms are used analogously to the synthetic setup and the resulting cumulative frequency of the position error over the grid agent positions is shown in Fig. 13. For comparison, we include our previous results from [32], where this evaluation was performed using an omni-directional antenna at the anchor. We identify a significant improvement, especially regarding estimation outliers (i.e. errors of more than half a meter). By consideration of the DMC viaĈ n (m) in the DMC-based method, the result can be improved even further. Specifically, the 80 % error is almost halved from 20 cm for the AWGN-based method down to 12 cm, and the amount of estimates to achieve an error of less than 10 cm goes up from 50 % to 72 %.
To gain further insight in the information provided by the SMCs and the beampatterns, we analyze |b m (φ k )| 2 SINR k , which is shown in Table 1. To obtain the SINRs, we apply (20)   Fig. 12) obtained using (20) with the estimated amplitudes from (48). The column added contains the sum over all antennas, whereas the entry omni contains the results from [32] using an omni-directional antenna in the same setup.
where we use the amplitude estimates obtained by (48). As seen in (38), these quantities indicate the quality of SMCs and highlight how each antenna can be used to focus on particular SMCs from different directions. E.g., we see that the SMC from the plasterboard east wall achieves a high SINR value from the east antenna measurement. For comparison, we include a sum of the weighted SINRs (similar to what is used in (38) to quantify the total information collected from each SMC with the multi-antenna system), as well as the SINRs of SMCs using a single omni-directional antenna at the anchor from [32] (the latter were obtained by methodof-moments estimation). The information gain of specific directional antennas is highlighted and the overall higher performance is justified. A variant of the AWGN-based algorithm was implemented on low-cost devices based on the DecaWave DWM1000 module as described in [31]. The system was evaluated in a field test and put in context with other indoor localization systems in terms of accuracy and required infrastructure.

VII. CONCLUSIONS
In this article, we investigated a multipath-based, singleanchor positioning system for indoor environments, which exploits non-coherent angular measurements from a set of directive antennas. Such measurements can be obtained for example from a mm-wave radio system with analog beamforming or a UWB transceiver with switched directive antennas. We derived and analyzed the position error bound (PEB), the Cramér-Rao lower bound (CRLB) on the position estimation error, for the proposed measurement system to gain insight about the achievable performance in comparison to a conventional antenna array requiring fully coherent processing of the antenna signals. Analysis of the PEB showed that the main contribution to high-accuracy (centimeter-level) positioning lies in the delays of specular multipath components (SMCs), whereas the directional antennas allow overlapping SMCs to be resolved, enhancing the robustness of the system. It was shown that non-coherent processing of directional measurements achieves a very similar performance compared to a fully coherent antenna array. Specifically, the bound on the RMS angle estimation error reduces by a factor of √ 3 when using the same number of antenna elements, whereas the delay estimation bound remains identical. This is a promising result, indicating that a high-accuracy, single-anchor positioning system exploiting angle information could be implemented efficiently with non-coherent (consecutive) measurements. Such measurements can be conducted with an adaptive, low-power analog antenna frontend, reducing the required number of radio chains to a minimum. In this regard, a position estimation algorithm was derived, showing the achievability of the PEB. Indeed, the performance is limited by dense multipath as selfinterference, not by measurement noise.

APPENDIX A DERIVATION OF EFIM GENERAL FIM
and (52) into (57), results in the EFIM for AoA estimation, given as where ξ angle k ∈ [0, 1] can be interpreted as the information loss on the SMC AoA estimation induced by the estimation of the complex amplitude α k in the DMC.

APPENDIX B ASSUMPTIONS
In the following, we will introduce assumptions for the SMC and DMC model that allow further insightful discussion in terms of contributions to position information.
• The DMC is uncorrelated between antennas, i.e., C v is block-diagonal (cf. (11)) with C (m,m ) v = 0 ∀m = m . For different antenna setups, there are different arguments supporting this assumption: For a conventional array consisting of omni-directional antennas, the DMC becomes approximately uncorrelated for an antenna-spacing of λ/2 and a uniform angular power spectrum [44]. In our considered case using directional antennas, it can be argued that each antenna at the anchor covers one sector in the azimuth plane with differently aligned main beam directions. As a result, the shapes of the beampatterns are approximately orthogonal such that With this, we may reduce the noise covariance matrix in (10) to block matrices on the main diagonal for which the mth matrix is given as