Member, IEEE Nil Garcia
email: nil.garcia@gmail.com

Member, IEEE Henk Wymeersch
email: henkw@chalmers.se

Fellow, IEEE Erik G Larsson
email: erik.g.larsson@liu.se

Fellow, IEEE Alexander M Haimovich
email: haimovic@njit.edu

Martial Coulon
email: martial.coulon@enseeiht.fr

Direct Localization for Massive MIMO

Keywords: MIMO, multipath channels, position measurement, 5G mobile communication, direction-of-arrival estimation, navigation, antenna arrays, signal processing algorithms, compressed sensing, sparse matrices, parameter estimation, base stations. I

Large-scale MIMO systems are well known for their advantages in communications, but they also have the potential for providing very accurate localization, thanks to their high angular resolution. A difficult problem arising indoors and outdoors is localizing users over multipath channels. Localization based on angle of arrival (AOA) generally involves a two-step procedure, where signals are first processed to obtain a user's AOA at different base stations, followed by triangulation to determine the user's position. In the presence of multipath, the performance of these methods is greatly degraded due to the inability to correctly detect and/or estimate the AOA of the line-of-sight (LOS) paths. To counter the limitations of this two-step procedure which is inherently suboptimal, we propose a direct localization approach in which the position of a user is localized by jointly processing the observations obtained at distributed massive MIMO base stations. Our approach is based on a novel compressed sensing framework that exploits channel properties to distinguish LOS from non-LOS signal paths, and leads to improved performance results compared to previous existing methods.

focused mainly on communications, it is also an enabler for high-accuracy localization [START_REF] Guerra | Position and orientation error bound for wideband massive antenna arrays[END_REF]. For instance, a finger-printing localization solution is proposed in [START_REF] Savic | Fingerprinting-based positioning in distributed massive MIMO systems[END_REF] for locating multiple users by means of distributed massive MIMO. A personal mobile radar with millimeter-wave massive arrays is proposed in [START_REF] Guidi | Personal mobile radars with millimeter-wave massive arrays for indoor mapping[END_REF] and used for simultaneous localization and mapping (SLAM) in [START_REF] Guidi | Millimeter-wave massive arrays for indoor SLAM[END_REF].

MIMO localization has received significant treatment in the technical literature, generally harnessing angle-of-arrival (AOA) estimation. Typically, a source emits a signal, and then in a two-step localization approach, the AOAs are measured at all base stations, and then, the source's location is found by triangulation. In benign open-air applications, where such methods are referred to as bearings-only target localization (BOTL) good performance can be observed [START_REF] Gavish | Performance analysis of bearing-only target location algorithms[END_REF]- [START_REF] Doganc | Passive emitter localization using weighted instrumental variables[END_REF]. However, in dense multipath environments, such as urban areas or inside buildings, the AOA estimates are biased in general. For that reason, pure AOA-based techniques [START_REF] Azzouzi | New measurement results for the localization of UHF RFID transponders using an angle of arrival (AOA) approach[END_REF] have not been very popular in harsh multipath environments, due to large localization errors [START_REF] Savic | Fingerprinting-based positioning in distributed massive MIMO systems[END_REF]. Massive arrays offer the possibility of precisely estimating the AOAs of the individual multipath components thanks to their high angular resolution [START_REF] Hu | An ESPRIT-based approach for 2-D localization of incoherently distributed sources in massive MIMO systems[END_REF], [START_REF] Lv | A beamspace approach for 2-D localization of incoherently distributed sources in massive MIMO systems[END_REF]. Nonetheless, measuring multiple AOAs at each base station requires identification of the AOA of the line-of-sight (LOS) paths. A possible strategy is to select the strongest arrival as LOS [START_REF] Klukas | Line-of-sight angle of arrival estimation in the outdoor multipath environment[END_REF]. However, the LOS path may be damped or obstructed as it is often the case indoors [START_REF] Spencer | Modeling the statistical time and angle of arrival characteristics of an indoor multipath channel[END_REF], [START_REF] Sen | Avoiding multipath to revive inbuilding wifi localization[END_REF]. Another option is to combine all estimated AOAs in a fusion center and perform data association [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF], but this is an NP-hard problem for which the optimal solution cannot be computed efficiently.

An alternative way to tackle localization problems, is to use a direct localization approach [START_REF] Weiss | Direct position determination of narrowband radio frequency transmitters[END_REF]. Contrary to traditional techniques, the location of the source is estimated directly from the data, without estimating intermediate parameters, such as the AOAs of the LOS paths. The concept of direct localization was first introduced in [START_REF] Wax | Location and the spectral density estimation of multiple sources[END_REF], [START_REF] Wax | Optimum localization of multiple sources by passive arrays[END_REF], and later applied to AOA-based localization [START_REF] Wax | Decentralized processing in sensor arrays[END_REF] and, more recently, to hybrid AOA-TOA (timeof-arrival) localization [START_REF] Weiss | Direct position determination of narrowband radio frequency transmitters[END_REF], [START_REF] Weiss | Direct position determination of multiple radio signals[END_REF]. However, all these methods were designed for pure LOS environments. Some direct localization techniques [START_REF] Chen | Maximum-likelihood source localization and unknown sensor location estimation for wideband signals in the near-field[END_REF], [START_REF] Bialer | Maximum-likelihood direct position estimation in dense multipath[END_REF] targeted to multipath scenarios exist in the literature, but they are not tailored to AOA information and massive arrays. A requirement of direct localization is that the signals, or a function of them, are sent to a fusion center that estimates the source's locations. In general, it is easier to achieve such a topology indoors as the distances are smaller. In the case of cellular networks, cloud radio access networks (C-RAN) [START_REF]C-RAN: The road towards green RAN[END_REF], [START_REF] Wu | Green Communications: Theoretical Fundamentals, Algorithms and Applications[END_REF] may provide the required infrastructure. C-RAN is a novel architecture for wireless cellular systems whereby the base stations relay the received signals to a central unit which performs all the baseband processing.

Our main contribution is a novel localization technique, called Direct Source Localization (DiSouL), that jointly processes the snapshots of data acquired at each base station in order to directly estimate the location of the source. Thanks to the high angular resolution of the massive arrays, the AOAs of the LOS paths can be used to precisely estimate the location of the source. DiSouL relies on a new compressive sensing framework 1 which exploits the fact that LOS components must originate from a common location whereas NLOS components have arbitrary AOAs. Contrary to previous AOA-based existing methods which assume that the LOS path is the strongest, this property enables DiSouL to infer the source position from the LOS AOAs even when the LOS paths are weaker than the NLOS paths. In comparison to ranging-based methods, DiSouL does not require large signal bandwidths in order to position a source accurately. To improve the signal-to-noise-ratio, DiS-ouL preprocesses the received signals with a matched filter, and then, samples at one time instant. We show how to determine such sampling instant at each base station based on a variation of the threshold matched filter [START_REF] Dardari | Threshold-based time-of-arrival estimators in UWB dense multipath channels[END_REF]. In addition, two mechanisms are presented to lower the computational burden and increase the accuracy. The first variation uses coarse TOA estimates at each base station to narrow the search area, while the second mechanism relies on a modified version of the grid refinement procedure [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. Finally, to validate the theory, numerous numerical results are provided, showing that DiSouL can achieve sub-meter accuracy with high probability which is sufficient for many applications (e.g., positioning of users in cellular networks, personal navigation, etcetera). In summary the contributions are: r A novel localization technique for massive MIMO sta- tions based on AOA information and assisted with TOA estimates.

r The mathematical formulation of a framework that enables the detection and estimation of LOS and NLOS paths without formulating a data-association problem.

r A simple generalization of the threshold matched filter for TOA estimation for arrays of antennas.

r A grid refinement procedure for lowering the computa- tional complexity of the proposed localization technique. Notation: • 1 , • 2 and • 2,1 denote the ℓ 1 -norm, ℓ 2norm and ℓ 2,1 -norm, respectively, and • 0 is the pseudo-ℓ 0norm which counts the number of non-zero elements.

II. SYSTEM MODEL

We consider a two-dimensional scenario with one user (or source) and L massive MIMO base stations equipped with arrays of S l antennas each. The user is located at p = [p x , p y ] T in an area R known a priori, the center of gravity of the stations' arrays are located at pl = [p x l , py l ] T and assumed to be in the far field with respect to the source. All arrays are equipped with fully digital processing, i.e., one radio frequency baseband chain per antenna [START_REF] Rusek | Scaling up MIMO: Opportunities and challenges with very large arrays[END_REF], [START_REF] Vieira | A flexible 100-antenna testbed for massive MIMO[END_REF]. We denote by a l (θ) the array response vector at base station l for a ray impinging with angle θ. The array response implicitly accounts for the array configuration and the antenna pattern of each antenna element.

The source broadcasts a known signal s(t) with half-power bandwidth B, which propagates through the multipath environment, resulting in a received signal at base station l given by

z l (t) = z LOS l (t) + z NLOS l (t) + n l (t) 0 ≤ t < T obs , (1)
where

z LOS l (t) = α l a l (θ l (p)) s(t -τ l (p)) (2)
z NLOS l (t) = P l m = 1 α m l a l (θ m l )s(t -τ m l ), (3)
in which T obs is the observation time, each component of n l (t) is white Gaussian noise with spectral density N 0 , α l is an unknown complex scalar, θ l and τ l are the angle of arrival (AOA) and time of arrival (TOA), all related to the line of sight (LOS) path, while α m l , θ m l , and τ m l are the channel gain, AOA, and TOA of the m-th NLOS component, for the P l NLOS path. All these parameters are unknown. The signal is narrowband with respect to the arrays, i.e., at each array l the amplitudes {α l }, {α m l } do not change across antennas. The LOS parameters τ l (p) and θ l (p) are related to the source position through

τ l (p) = p -pl /c (4) θ l (p) = arctan p y -py l p x -px l + π • 1 (p x < px l ) , ( 5
)
where c is the speed of light, while the range of the arctangent function is -π/2 ≤ arctan(x) < +π/2, the angle is computed with respect to the x-axis and anticlockwise, and 1(P) is one if the logical expression P is true 2 . Let pl,a = [p x l,a , py l,a ] T be the position of antenna a at BS l, relative to the array's center of gravity. For arrays without mutual antenna coupling and isotropic antennas, the array response a l (θ) for a given AOA θ admits the following close-form expression

[a l (θ)] a = exp 2πi λ pT l,a cos(θ) sin(θ) , (6)
where [•] a denotes the a-th component, i is the imaginary unit and λ is the wavelength of the carrier. In practice, for non-ideal arrays with mutual coupling and different antenna gains, a l (θ) is not computed mathematically but is measured during the array calibration process. The collection of vectors a l (θ) for different directions θ is often referred to as the array manifold. Define the emitted signal autocorrelation r s (t) = T obs 0 s * (τt)s(τ ) dτ . We generate a discrete-time obser- 2 The summand π • 1 p x < px l is added for resolving the ambiguity caused by the fact that arctan(y/x) = arctan(-y/-x). α m l r s (tτ m l ) a l (θ m l ) + n MF l (t), [START_REF] Guidi | Millimeter-wave massive arrays for indoor SLAM[END_REF] and sampling at time t l ,

zl = z MF l (t l ) = ᾱl a l (θ l (p)) + P l m = 1 ᾱm l a l (θ m l ) + nl (8)
where ᾱl = r s (t lτ l (p))α l , ᾱm l = r s (t lτ m l )α m l and nl ∼ CN (0, σ 2 I) where σ 2 = N 0 because, without loss of generalization, the pulse energy is normalized to one (r s (0) = |s(τ )| 2 dτ = 1). The signals in [START_REF] Gavish | Performance analysis of bearing-only target location algorithms[END_REF] are the input of the proposed method. Contrary to beamforming applications, and similar to direction-of-arrival techniques, the observations across each array zl are not weighted and linearly combined into one output. Such weights would fix the array beampattern and are usually designed to make the array point into one or multiple directions. In our case, since no beamforming is performed, the arrays do not have a "favored direction". In regards to the source, we assume that it has an omnidirectional antenna, in which case the transmitted energy would be the same towards all BSs.

In order to ensure that ᾱl = 0, the signals must be sampled at a time where the energy of the LOS pulse is not zero (i.e. while r s (t lτ l (p)) = 0). In addition to the sampling times, we will also compute an upper bound on the TOA of the LOS paths at each base station that will enhance the proposed method, which for brevity are simply called TOA estimates. In a nutshell, the objective of our work is to determine the sampling times and the TOA estimates from {z MF l (t)} L l= 1 , and then, determine p from {z l } L l= 1 .

III. PROPOSED METHOD

A. Principle

The proposed method exploits the high angular resolution of massive arrays, enabling detection and estimation of the AOAs of the distinct multipath arrivals. The TOA estimates are not used for precise localization, but rather as constraints, limiting the source location in a convex set (see Fig. 1 for a visual example). By processing the snapshots zl at all base stations jointly, we are able to separate the LOS paths from the NLOS paths. The proposed method exploits the information carried by the LOS paths in order to position the user with high accuracy. The implicit assumption here is that, in each BS, the LOS path is present. A method for coping with scenarios where some BS are in NLOS is explained in Section IV-C. Roughly speaking, the procedure of DiSouL is as follows. Determine a coarse positively biased estimate of the TOA at each base station, and use these to narrow the search area to a convex set containing the source. Then, using the signal model ( 8), we formulate a convex optimization problem, directly providing an estimate of p. In contrast to indirect approaches, where AOAs are estimated first and the source location is determined afterwards, we do not have to deal with an NP-hard data-association problem. We bypass this problem, as only position p is estimated and not AOAs of the LOS paths. While the TOA estimates need to be an upper bound on LOS TOAs, i.e. positively biased, the sampling should happen at an instant where the energy of the LOS arrival is maximized with respect to the energy from the NLOS arrivals. Thus, in general the sampling times will be smaller than the TOA estimates. The proposed technique is designed for a single user, although it could possibly be extended to multiple users. In systems operating under an orthogonal multiple access (OMA) scheme (e.g., TDMA, FDMA, synchronous CDMA), the signals of each user could be filtered out from the rest and treated separately.

B. TOA Assistance

A delay estimation technique is described for computing coarse TOA estimates {τ l } L l= 1 , which are then used to reduce the search area. Nonetheless, other delay estimation techniques are possible and may lead to better results. We rely on a generalization of the threshold MF [START_REF] Dardari | Threshold-based time-of-arrival estimators in UWB dense multipath channels[END_REF] for multiple antennas. This technique has the advantage of being relatively simple, but more importantly, in the presence of multipath, due to the fact that multipath components overlap in time, it generates positively biased TOA estimates as required by our method. Note that the poor time resolution is not a concern because the goal of the proposed method is to resolve the LOS from the NLOS components in the 'angular domain' rather than the time domain. The threshold MF first applies a matched filter to the received signals [START_REF] Guidi | Millimeter-wave massive arrays for indoor SLAM[END_REF], and then seeks the first peak that exceeds a properly chosen threshold. By ignoring the part of the signals below the threshold, the threshold MF avoids false detections due to the noise, and by selecting the first peak above the threshold, it obtains a better estimate of the LOS TOA than, for instance, selecting the strongest peak. We propose a simple generalization of this approach from the single-antenna case to the case of arrays with multiple antennas.

Let z NC l (τ ) be the non-coherent aggregation of the observed signals at all antennas after matched filtering:

z NC l (t) = z MF l (t) 2 2 , (9)
where z MF l (t) was introduced in ( 7). The TOAs {τ l } L l= 1 are estimated by selecting the first peak4 that exceeds a certain threshold:

τl = find-1st-peak z NC l (t) : z NC l (t) ≥ η . (10)
In practice, values of z NC l (t) may only be available at discrete instants. In such case the location of the peak may be obtained by parabolic fitting [START_REF] Cespedes | Methods for estimation of subsample time delays of digitized echo signals[END_REF]. Following [START_REF] Dardari | Threshold-based time-of-arrival estimators in UWB dense multipath channels[END_REF], the threshold is selected so that the probability of early false alarm is very low. An early false alarm event is defined as detecting a peak due to noise before the true TOA of the LOS path. Since the threshold MF ignores second and later peaks (above the threshold), the following expression on the probability of false alarms only considers those peaks happening before the LOS TOA. When the false alarm probability is small, it can be approximated as [START_REF] Dardari | Threshold-based time-of-arrival estimators in UWB dense multipath channels[END_REF]

P FA ≈ 1 + (1 -q noise ) N TOA -1 N TOA q noise , (11)
where N TOA = T obs /T corr and T corr is the waveform correlation time, which for most types of waveforms, is well approximated by the inverse of the bandwidth T corr = 1/B. The term q noise is defined as the probability that the noise in absence of signal exceeds threshold η

q noise = Prob n MF l (t) 2 2 ≥ η = 1 -F 2S l 2η σ 2 , ( 12
)
where S l is the number of antennas at BS l, σ 2 is the noise variance at each antenna, and F k (x) be the cumulative distribution function of the Chi-squared distribution with k degrees of freedom evaluated at x. The threshold η resulting in the desired P FA can be found by performing a one-dimensional search of [START_REF] Azzouzi | New measurement results for the localization of UHF RFID transponders using an angle of arrival (AOA) approach[END_REF]. Assuming all TOA estimates are positively biased, then we can create a set

F = {π ∈ R 2 : π -p l 2 ≤ c τl , ∀l} (13)
and narrow the search area to R ∩ F instead of R. In the unlikely event that not all TOA estimates are positively biased, it is possible that F = ∅. In such a case, we v expand F = by increasing all TOA estimates by a constant value until F ∅. We have chosen v = 1/B, where B is the signal bandwidth, though the value of v turns out to not be critical for the localization performance. Reducing the search area has two benefits. On one hand it lowers the computational complexity of the optimization problem that is proposed in the next section. On the other hand it can positively affect the accuracy of the localization procedure.

C. Sampling Times

The outputs of the matched filters ( 8) at all base stations ( 8) are sampled at times {t l } L l= 1 . Contrary to the estimation of the TOAs {τ l } L l= 1 , the goal is to sample {z MF l (t)} L l= 1 at an instant where there is as little as possible NLOS interference and as much as possible energy from the LOS component. One solution would be to sample at the same time than the TOA estimates, however, due to their bias, we will rarely sample at the time of maximum LOS energy, and moreover, the NLOS interference may be very large as well. We propose an alternative strategy, which has been verified numerically to lead to better results in Section VI for the propagation conditions under consideration. The idea is to use the same threshold matched filter for TOA estimation, but instead of selecting the time of the first peak, we select the instant when the received signal crosses the threshold for the first time, i.e.,

t l = min t : z NC l (t) ≥ η . ( 14
)
Comparing the expression of the sampling time to that of the TOA estimate [START_REF] Doganc | Passive emitter localization using weighted instrumental variables[END_REF], it follows that t l ≤ τl for all l.

D. Localization

To solve the localization problem, we rely on tools from compressive sensing. Specifically, we propose a grid-based approximate solution to the problem of localizing a source on a continuous map which exploits the notion of sparsity and rowsparsity [START_REF] Tropp | Algorithms for simultaneous sparse approximation[END_REF]. To this end, first we introduce a uniform grid of Q locations

L = {π 1 , . . . , π Q } ⊂ R ∩ F, (15)
and a uniform grid of M l angles for each base station array

A l = {ϑ 1 , . . . , ϑ M l } ⊂ [0, 2π). ( 16
)
The main assumption here is that the source is positioned on a grid location, and that the AOAs of the NLOS paths are also in the grid of angles. While beamforming suffers from the Rayleigh resolution limit [START_REF] Smith | Statistical resolution limits and the complexified Cramér-Rao bound[END_REF], which is independent of the SNR, an advantage of sparsity-based techniques is that they can resolve multipath components within a Rayleigh cell (i.e., achieve super-resolution) [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF], [START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF], [START_REF] Fuchs | Multipath time-delay detection and estimation[END_REF]. Let X ∈ C Q ×L be a matrix whose entry on row q, column l is denoted by x q l and represents the complex gain of a LOS path from grid location π q to base station l. Let y m l be the complex gain of a NLOS path arriving at the l-th base station with angle ϑ m . Then, by definition, only one row in X is different from zero, and y m l = 0 only if ϑ m is equal to the AOA of a NLOS path at base station l. Thus, if the grids are dense enough, X is rowsparse and y l is sparse for all l. It is well known in the compressive sensing literature [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF], that row sparsity can be induced by

minimizing the ℓ 2,1 -norm, i.e. X 2,1 = Q q = 1 L l= 1 |x q l | 2 ,
and that sparsity can be induced by minimizing the ℓ 1 -norm, i.e. y l 1 = M l m = 1 |y m l | where y l = [y 1l , . . . , y M l l ] T . Thus, with all this in mind, we propose to solve the following optimization problem

min X,y l w X 2,1 + L l= 1 y l 1 (17a) s.t. L l= 1 zl -ẑl 2 2 ≤ ǫ (17b) ẑl = Q q = 1 x q l a l (θ l (π q )) + M m = 1 y m l a l (ϑ m ) , ∀l. (17c)
This is a second-order cone program (SOCP) for which very efficient solvers exist. The optimization variables are X and {y l } L l= 1 . The vector ẑl as defined in (17c) is a reconstruction of zl for a given choice of X and {y l } L l= 1 . The parameter ǫ establishes the maximum allowed mismatch between the observations and the reconstruction.

Intuitively, problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] seeks the sparsest number of source locations and NLOS paths that can describe the observations {z l }. Precisely, minimizing X 2,1 , induces a sparse number of source locations (hopefully a single location), and minimizing L l= 1 y l 1 , induces a sparse number of NLOS paths. These are two different types of sparsity, and therefore, in order to combine them together some kind of balancing coefficient or weight w needs to be added. Intuitively, the parameter w ensures that not all signal energy is assigned to only LOS or NLOS components. In particular, if w = 0, the objective function is not penalized by filling matrix X with non-zeros, thus, the reconstruction will always favor a representation based on LOS paths rather than NLOS paths. On the contrary, if w → ∞, then the objective function will highly penalize any non-zero in matrix X, thus, the reconstruction will always favor a representation based on NLOS paths. Thus, it is clear than in our multipath scenario, consisting of LOS and NLOS paths, if a suitable w exists, it must be larger than 0 and finite. Below suitable choices for ǫ and w are proposed.

Remark: While this technique searches a source on a plane, it can be generalized to a three-dimensional search, at a cost of computational complexity. It is also possible that the technique may improve its robustness against multipath because in two dimensions two distinct NLOS bearing lines always intersect but in three dimensions they generally do not.

IV. PARAMETER SELECTION

In this section, we motivate the choice for ǫ and w in the optimization problem of Section III-D. The choices do not guarantee recovery of the correct position and are derived under simplified assumptions.

A. Setting the Parameter ǫ

The parameter ǫ in [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] defines the allowed mismatch between the observations and the reconstruction. Typically, ǫ is a bound on the noise. Since the noise is Gaussian, it is unbounded, and instead ǫ is chosen so that it is a bound with high probability, i.e.

Prob L l= 1 nl 2 ≤ ǫ = γ. ( 18
)
where γ is, for example, 0.99. Because nl are random white Gaussian vectors of length S l , it follows that the error normalized by the noise variance 2σ -2 L l= 1 nl 2 is a Chi-squared random variable with 2 L l= 1 S l degrees of freedom. Let F k (x) be the cumulative distribution function of the Chi-squared distribution with k degrees of freedom evaluated at x and F -1 k (y) its inverse function evaluated at y. Then, the value of ǫ can be computed as

ǫ = σ 2 2 F -1 2 L l = 1 S l (γ) . (19)
In low SNR conditions, it is possible that the aggregated energy of all snapshots is not larger than ǫ, i.e.

L l= 1 zl 2 2 ≤ ǫ, (20)
making problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] have the trivial all-zeros solution, thus, failing to estimate the location of the source. In this case we propose to simply look up the location whose LOS components correlate the most with the snapshots:

p = arg max π∈L L l= 1 a H l (θ l (π)) zl 2 a l (θ l (π)) 2 2 . ( 21
)

B. Setting the Parameter w

In order to obtain an expression for w, we will not prove that the AOAs are correctly recovered by solving [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF], but rather that, under proper selection of w, if the AOAs are correctly recovered, then they can also be correctly identified as either LOS or NLOS. The key property that will dictate the value of w, and in turn estimate the correct source location is based on the following definition.

Definition 1 (consistent location): A location π is consistent with L paths, if the AOAs of the direct paths between π and the base stations are true AOAs, i.e.

θ l (π) ∈ Θ l for l = 1, . . . , L, (22)
where Θ l is the set of true AOAs at base station l

Θ l = θ l (p), θ 1 l , . . . , θ P l l . ( 23
)
By definition, the true source location p is consistent with the L paths because the LOS components travel in a straight line. To find a criterion for the weight, we restrict ourselves to a simplified version of the problem and then later evaluate the criterion in a more realistic setting. Our analysis on the weight criterion is limited through the three following assumptions. A1) Besides the source location p, no other location is consistent with L paths. A2) The grids L and {A l } L l= 1 are sufficiently dense to contain the source location p and all AOAs, respectively. A3) Denoting by Θl the estimated AOAs at base station l, i.e.

Θl = {θ l (π q ) : x q l = 0} ∪ {ϑ m l : y m l = 0} , (

then Θl = Θ l , ∀l. In other words, the solution of ( 17) is able to recover the true AOAs. This assumption is reasonable in high SNR conditions. These assumptions are only used for deriving a suitable value of w. In practice, in realistic multipath environments, some of these assumptions may break down. Thus, in Section VI, the proposed method is tested versus a realistic indoor multipath channel model, and shown to recover the user position with sub-meter accuracy and high probability.

Lemma 1: Assume A2) and A3). If w > √ L -1, then any estimated location output by problem ( 17) is consistent with L paths (in the sense of Definition 1).

Proof: See Appendix A. Lemma 2: Assume A2) and A3). If w < √ L, then problem (17) outputs at least one location (i.e. X = 0).

Proof: See Appendix B.

The two lemmas lead directly to the following theorem, which guarantees the correct recovery of the source location.

Theorem 1: If Assumptions A1), A2), and A3) hold, then a sufficient condition for the correct recovery of the source location is

√ L -1 < w < √ L. ( 25
)
Proof: If w < √ L, by Lemma 2 al least one estimated location is output by problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF]. Moreover, if w > √ L -1, by Lemma 1 any estimated location is consistent with L paths. However, according to Assumption A1), only the location of the source is consistent with L paths, thus completing the proof.

C. The Case of Non-LOS

In practice, BSs may be in non-line-of-sight (NLOS) because the LOS paths are blocked or attenuated (making the LOS path pass undetected). The proposed technique relies on the presence of the L LOS paths for achieving high precision localization, and it may break down when the base stations are in NLOS. If the weight in (17a) is chosen according to Lemma 1, any location estimate output by problem (17) must be consistent with L paths. However, if one base station is in NLOS, then the source will only be consistent with L -1 paths, and therefore, the location of the source will not be a solution to [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF]. Thus, adjusting the weight requires a priori knowledge of the number of LOS base stations. We can adjust the weight as follows. Let L * be the number of base stations in LOS with the source, and let L be an estimate of L * . Obviously, L * ≤ L. Furthermore, assume no other location besides the location of the source is consistent with L * paths. We start by assuming that all base stations are in LOS, i.e. L = L, set the weight according to solve [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] to obtain X and y l , ∀l 7:

if (X ≡ 0) then 8:

L ← L -1 9: else 10: q = arg max q X q ,: 2 11: p = π q 12: end if 13:

end while 14: else 15: estimate p by ( 21) 16: end if Theorem 1, and solve problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF]. According to Lemma 1, the location of the source will be estimated only if it is consistent with L paths. If L * < L, the solver will return X = 0. When this event is detected, L can be reduced and ( 17) solved again. This procedure can be repeated as shown in Algorithm 1. Note that if X = 0, the location with strongest gains is returned as the estimate p (see lines [START_REF] Doganc | Passive emitter localization using weighted instrumental variables[END_REF][START_REF] Azzouzi | New measurement results for the localization of UHF RFID transponders using an angle of arrival (AOA) approach[END_REF]. Theoretically speaking, since two non-parallel paths always intersect on the plane, to ensure that we can distinguish LOS from NLOS paths, and in turn, recover the user position, the number of LOS BSs should be L * ≥ 3. However, when the LOS paths are the strongest, we have verified numerically, that the user can be recovered with L * = 2 as well.

Based on the choices for ǫ and w, Algorithm 1 summarizes the proposed solution strategy.

V. GRID REFINEMENT

Dense grids of locations and angles are necessary to achieve fine resolution, but making the grids too dense results in large computation time. The computational complexity of solving [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] scales as O((QL + l M l ) 3.5 ) [START_REF] Lobo | Applications of second-order cone programming[END_REF], where Q and M are the number grid locations and angles respectively, and L is the number of base stations. This motivates an adaptive gridrefinement strategy originally proposed in [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. The idea behind the grid refinement approach is to start with a coarse grid of locations and angles; subsequently, the grid is refined around the estimated locations and angles and the optimization problem ( 17) is solved again. This procedure can be repeated until a certain grid resolution has been achieved or a stopping criterion has been met. Thus, the benefits of grid refinement are two-fold: lower computational complexity and fine grid resolution.

In comparison to previous grid refinement approaches [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF], [START_REF] Hyder | Direction-of-arrival estimation using a mixed norm approximation[END_REF], ours is more complex due to the two different types of grids used to describe the observed data. At iteration k of the grid refinement process, we will denote the position grid by L (k ) and the angle grid (for base station l) by A (k ) l . At iteration k = 0, the grids are uniform over R ∩ F and [0, 2π), respectively. The resolutions in L (0) and A l (0) are set to π res ∈ R + and ϑ res ∈ R + , respectively. Define the following operators:

grid( Π, δ) = {π ∈ R 2 : π = π + [i˜j] T δ, π ∈ Π, i, j ∈ {-2, -1, 0, 1, 2}} (26)
grid( Θl , δ) = {θ ∈ [0, 2π) : θ = θ + iδ, θ ∈ Θl , i ∈ {-2, -1, 0, 1, 2}}. (27)
We can then set L (k ) = grid( Π, π res /2 k ) and

A (k ) l = grid( Θl , ϑ res /2 k ) ∪ {[θ l (π)] ϑ res /2 k : π ∈ Π}
, where [x] y rounds x to the nearest multiple of y. Each successive grid of locations and angles includes the estimated points and their neighboring points. In this case we have chosen to include twenty-four and four neighbor points for the position grid and angle grid, respectively, but other choices of neighbors are possible as well. In addition, the grid of angles also incorporates the angles related to the estimated locations. It has been empirically verified that this is necessary for the correct performance of this grid refinement approach.

Because at each step the previously estimated points are included in the next grid, the solution at step k is a feasible solution at step k + 1. This ensures that the optimum value of the optimization problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] cannot increase as iterations progress. Since the objective function is bounded from below by zero, by the monotone convergence theorem [START_REF] Schechter | Handbook of Analysis and Its Foundations[END_REF], the grid refinement procedure must converge. In practice, the refinement process is halted when the progress between two consecutive steps is negligible. Denote as f (k ) opt the optimum value of problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] at step k, then the grid refinement is stopped at step k if

f (k -1) opt -f (k ) opt f (k -1) opt < β, ( 28
)
where β is a small value, e.g., β = 10 -3 .

A. The DiSouL Algorithm

The summary of the DiSouL algorithm is now presented, comprising the basic Algorithm 1, as well as the TOA assistance of Section III-B and the grid refinement of Section V.

VI. NUMERICAL RESULTS

In this section, we illustrate the performance of the localization method and compare it to other existing techniques. Unless otherwise stated, all numerical examples are run using the following parameters. The source is positioned randomly within an area of size 100 × 100 m. Four base stations are positioned at the corners; if the origin of the coordinate system is taken to be in the middle of the area, the base stations are at coordinates [45 m, 45 [START_REF] Fuchs | On the application of the global matched filter to DOA estimation with uniform circular arrays[END_REF] with a radius of 24 cm because it makes the inter antenna spacing equal to half wavelength. We opt for UCAs instead of uniform linear arrays (ULAs) because for the same number of antennas and equal inter antenna spacing, their far field region [START_REF]Standard Test Procedures for Antennas[END_REF] starts at a much shorter distance (11 m for a UCA and 51 m for a ULA), and contrary to ULAs, UCAs have equal angular resolution towards all directions on the plane. Also, we simulate four base stations because it makes positioning a bit more robust to multipath than using only three base stations 5 while keeping the number of BSs to a reasonable number. The initial grid resolutions of DiSouL are π res = 5 m and ϑ res = 5.71 • .

A. Validation of Theorem 1

To illustrate Theorem 1, we synthesize a set of snapshots [START_REF] Gavish | Performance analysis of bearing-only target location algorithms[END_REF] according to the scenario plotted in Fig. 2 and ignore any time delay information. The source is positioned at [18 m,31 m] and a reflector is positioned at [25 m,-7 m]. As visualized in Fig. 2, all base stations receive a LOS component, and except for the top left base station, they also receive a NLOS component bounced from the reflector. From Fig. 2, it is apparent that the source location is consistent with 4 paths, the reflector is consistent with 3 paths, and all other locations in that area are consistent with 2 paths or less. We hypothesize that for a sufficiently fine grid, and for a sufficiently high SNR, the probability of recovering the correct source location will be high if the weight is picked according to Theorem 1. The SNR l for the snapshots ( 8) is defined as

SNR l = S l E{|ᾱ l | 2 }/E{ nl 2 2 } = E{|ᾱ l | 2 }/σ 2
and is equal for all base stations SNR l = SNR. Every data point in Fig. 3 is generated by running 100 Monte Carlo runs, where at each run the signal strengths and phases of all multipath components are randomized according to Rayleigh and uniform distributions, respectively. The location of the source is estimated by running Algorithm 2, wherein the solution to optimization problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] is obtained by the solver Mosek [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual, version 7.1[END_REF]. Fig. 3 plots the empirical probability that the localization error is smaller than 1m as a function of w 2 . According to Theorem 1, a sufficient condition for recovering the location of the source is that the square of the weight satisfies L -1 ≤ w 2 ≤ L. The figure shows that, in this case, for L = 4, the range of values Algorithm 2: Grid refinement.

1: given a coarse grids of locations L (0) and angles A

l , ∀l 2: set k = 0 3: while (28) not satisfied do 4:

solve [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] with

L = L (k ) and A l = A (k ) l 5: extract locations Π = {π (k ) q ∈ L (k ) : x (k ) q 2 = 0} 6: extract angles Θl = {ϑ (k ) m l ∈ A (k ) l : y (k )
m l = 0}, ∀l 7: increase k 8:

set

L (k ) = grid( Π, π res /2 k ) 9:
trim grid of locations L (k ) through TOA assistance 10: set

A (k ) l = grid( Θl , ϑ res /2 k ) ∪ {[θ l (π)] ϑ res /2 k : π ∈ Π}
11: end while Algorithm 3: DiSouL.

1: set η using [START_REF] Azzouzi | New measurement results for the localization of UHF RFID transponders using an angle of arrival (AOA) approach[END_REF] for the desired P FA (e.g., P FA = 10 -2 ) 2: estimate TOAs {τ l } L l= 1 using (10) 3: create initial grid of locations L (0) and angles A (0) l , ∀l 4: trim grid of locations through TOA assistance 5: compute sampling times {t l } L l= 1 using 6: obtain the snapshots of data by applying MF and sampling at instants {t l } L l= 1 as in (8) 7: estimate source location p by Algorithm 1 where line 1 is replaced with Algorithm 2 [START_REF] Chin | Emerging technologies and research challenges for 5G wireless networks[END_REF][START_REF] Guerra | Position and orientation error bound for wideband massive antenna arrays[END_REF] yields the correct source location with sub-meter accuracy with probability 1 for a sufficient high SNR.

w 2 ∈ [

B. Localization Performance in Realistic Multipath Channel

In this section, we perform Monte Carlo simulations in which DiSouL is compared to indirect and direct localization techniques:

1) SR-LS [START_REF] Beck | Exact and approximate solutions of source localization problems[END_REF], using TOA obtained by the time delay estimator of Section III-B. SR-LS takes as inputs a set of ranges (related to the TOAs by the speed of light) obtained at distributed BSs and outputs the user position by solving a least-squares (LS) fit. Because the maximum likelihood (ML) solution does not lead to a convex optimization problem, the LS solution is only an approximation. 2) IV [START_REF] Doganc | Passive emitter localization using weighted instrumental variables[END_REF], using AOA information, obtained by applying beamforming [START_REF] Godara | Application of antenna arrays to mobile communications. II. Beam-forming and direction-of-arrival considerations[END_REF] on the snapshots (8) and selecting the angle associated with the strongest peak. IV is a closedform estimator that uses the AOAs measured at distributed BSs to triangulate the user position. In pure LOS environments, IV is consistent and asymptotically (with the number of BSs) unbiased.

3) The Stansfield estimator [START_REF] Gavish | Performance analysis of bearing-only target location algorithms[END_REF], using hybrid TOA-AOA. The Stansfield estimator is sometimes described as an AOAbased only estimator. It uses coarse range estimates in order to approximate the AOA-based ML estimator (in absence of multipath) by a closed-form solution. Its formula is a refined version of IV. 4) DPD [START_REF] Weiss | Direct position determination of narrowband radio frequency transmitters[END_REF], a direct localization hybrid TOA-AOA technique, operating directly on the received signals ( 1)-( 3). DPD is essentially the ML estimator for a pure LOS environment when directly operating on the received signals. Since the received signals depend on the user position through their AOAs and TOAs, DPD is the optimal hybrid technique in absence of multipath. In practice it requires computing a figure of merit for each grid location on the map, and the one yielding the largest value is the user's position estimate. Due to this two-dimensional search, it is, in general, more computationally expensive than indirect techniques. The AOAs for IV and the Stansfield estimator are estimated by performing beamforming. More sophisticated AOA estimation techniques, such as MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], are not applicable because they require multiple snapshots and break down in the presence of multiple correlated arrivals such as is the case of multipath. A high precision alternative to beamforming is ℓ 1 -SVD [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. However, we have observed in our numerical results that ℓ 1 -SVD performs similar to beamforming due to the fact that the AOA estimation errors are caused by peak ambiguities and not the lack of angular resolution. Thus, errors happen mostly when the LOS component is attenuated or blocked by obstacles.

The source emits a Gaussian pulse s(t) at 7 GHz carrier frequency. We simulate the received signal at each antenna after down-conversion to baseband and sampling. An oversampling factor of 3 is used. It is a assumed a half power bandwidth of B = 30 MHz and 10 log 10 E N 0 = 10 dB, where E = E |α l | 2 is the energy of the received LOS component before sampling (2) (same energy for all l) and N 0 is the noise spectral density at each antenna.

Every data point in the figures is generated by running 1000 Monte Carlo runs. The parameters randomized at each Monte Carlo run are the multipath channels, the noise and the position of the user. The parameters of the multipath channel (1) are generated according to the statistical indoor channel model in [START_REF] Spencer | Modeling the statistical time and angle of arrival characteristics of an indoor multipath channel[END_REF], the user position is drawn from the uniform distribution over the search area, and the noise is independent and identically Gaussian distributed. In regards to the indoor multipath channel, the parameters' values are those of the Clyde building: cluster decay rate is 34ns, ray decay rate is 29ns, cluster arrival rate is 1/17ns, ray arrival rate is 1/5ns and angular variance is 26 • . On the average, at every base station, 99.9% of the energy in the snapshot [START_REF] Gavish | Performance analysis of bearing-only target location algorithms[END_REF] is contained in 8 discrete multipath arrivals, which, in general, have closely spaced AOAs.

Fig. 4 plots the cumulative density function of the localization error. Clearly, DiSouL achieves high precision accuracy with high probability, followed by DPD and the two-step approaches. To gain a more in-depth understanding, we will focus on the performance of the estimators at sub-meter errors, as a function of E/N 0 , bandwidth, number of antennas, channel properties, and calibration errors.

In Fig. 5, the probability of sub-meter precision is shown as a function of E/N 0 . Note that DiSouL outperforms all other techniques for most E/N 0 values. The TOA-based SR-LS performs poorly due to the positive bias of the TOA estimates. The AOA-based estimators can slightly improve on this performance, but are still worse than both direct localization approaches. As E/N 0 increases, we sample the snapshots at the time of crossing a threshold rather than at the peak (see Section III-B), which reduces the amount of NLOS multipath components that are included into the snapshots, but the resulting ratio between LOS energy and noise is more or less independent of E/N 0 . Thus, the benefit of increased E/N 0 is that we detect the signals sooner, thus diminishing the number of NLOS components making into the snapshots. In the low SNR regime (less than -5 dB in the figure), DiSouL is outperformed by DPD because the threshold matched filter used for sampling the signals and estimating the TOA may fail to detect any signal.

Fig. 6 plots the probability of sub-meter accuracy versus signal bandwidth. All techniques benefit from an increase of bandwidth. On the one hand, it is well known that a larger bandwidth results in better TOA estimates. On the other hand, since the pulse width is inversely proportional to the bandwidth, a larger bandwidth results in a shorter pulse. Hence, fewer NLOS multipath components are included into the snapshots (8), thus, decreasing the risk of errors in the AOA estimation. Fig. 7 evaluates the probability of sub-meter accuracy versus the number of antennas in each base station. To keep the separation between two neighboring antennas in each UCA equal to λ 2 , the diameter of the array in BS l is increased with the number of antennas according to the formula λ 2 sin(π /S l ) . For ten or more antennas, the expression of the diameter is well approximated by λS l 2π , meaning that the array size grows proportionally to the number of antennas S l , which is directly related to the angular resolution. This improvement in angular resolution allows DiSouL to resolve more multipath arrivals more precisely, and consequently, improve its localization accuracy as observed in Fig. 7. Between 10 and 80 antennas, DiSouL's probability of sub-meter accuracy improves linearly with the number of antennas, and then it saturates. A perfect probability of 1 is not achieved because the assumptions made by DiSouL (see Section IV-B) may not always hold true in the simulated multipath channel. On the contrary, the probability of sub-meter accuracy for the indirect techniques remains approximately the same. In particular, because SR-LS is purely TOA-based, the improvement in angular resolution has no impact. The other two indirect techniques, IV and Stansfield, improve very little because most of their errors are due to selection of the wrong path as LOS.

In Fig. 8, we tune some of the channel parameters controlling the rate of arrivals. In the statistical multipath channel model of [START_REF] Spencer | Modeling the statistical time and angle of arrival characteristics of an indoor multipath channel[END_REF], the times of arrival of the NLOS components are modeled by two parameters: the cluster arrival rate λ and the ray arrival rate λ. The measured values for the Clyde building of these two parameters were 1/λ = 17 ns and 1/λ = 5 ns. In order to study the localization accuracy as a function of the ray arrival time, in Fig. 8, 1/λ is varied between 5 ps and 5 µs while λ = 5 17 λ. As the ray inter-arrival time increases, the multipath channel becomes less dense. For very high inter-arrival times, the channel can be considered almost pure LOS, and as expected all techniques improve their localization accuracy.

In Fig. 9, we study the effect of calibration errors on the arrays, more specifically, the degradation in localization accuracy due to mismatches between the true and nominal gains/phases of the antennas. Based on the model in [START_REF] Friedlander | A sensitivity analysis of the MUSIC algorithm[END_REF], the true array response ãl is given by

[ã l (θ)] a = g l,a e iφ l , a [a l (θ)] a , (29)
where a l is the nominal (i.e., the array response used by DiSouL) array response [START_REF] Guidi | Personal mobile radars with millimeter-wave massive arrays for indoor mapping[END_REF], and g l,a ≥ 0 and φ l,a represent the unknown gains and phases, respectively, of the antennas. The gains (in

TABLE I AVERAGE EXECUTION TIMES OF THE MULTIPLE LOCALIZATION TECHNIQUES

dB) and phases of all antennas are drawn independently from Gaussian distributions, i.e. 20 log 10 g l,a ∼ N (ḡ dB , σ 2 dB ) and φ l,a ∼ N (0, σ 2 φ ), where ḡdB = -ln(10) 40 σ 2 dB because it makes E g l,a = 1. When the nominal array response matches the true one (i.e., σ dB = σ φ = 0), the probability of sub-meter accuracy is approximately 0.8. Obviously, as larger gain and phase errors occur at all antennas, the performance of DiSouL degrades. The probability of sub-meter accuracy drops from 0.8 to 0.7 at 38 • phase standard deviation (in the absence of gain mismatches), or at 6 dB gain standard deviation (in the absence of phase mismatches). Thus, DiSouL is remarkably robust to calibration errors.

Lastly, Table I plots the execution times of all techniques in a regular 3.6 GHz desktop computer. All stations are equipped with 70-antenna arrays. Due to the joint processing of the data at all base stations, the execution times of the direct techniques are much larger than those of the indirect techniques. In particular, DiSouL is substantially more computationally intensive than the other techniques because it needs to solve a relatively large optimization problem multiple times. Instead of using an off-the-shelf solver [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual, version 7.1[END_REF], it may be worth developing/using an algorithm that exploits the sparsity of the signals.

VII. CONCLUSIONS

This paper tackled the problem of narrowband localization in the presence of multipath, through a direct localization approach in a massive MIMO setting. We propose an original compressive sensing approach for the localization of sources emitting known narrow-band signals. Due to the high angular resolution of massive arrays, it is possible to estimate the AOAs of the multipath components. By jointly processing snapshots of several widely distributed arrays, we are able to estimate the source location precisely without explicitly estimating the LOS AOAs, and therefore, avoiding the challenging data association problem. The proposed technique, called DiSouL, achieves sub-meter localization with high probability in dense multipath environments with narrow-band signals. DiSouL requires no statistical channel knowledge except for the noise variance and therefore it is suitable for any multipath environment. Coarse TOA estimates at each array are used to reduce the execution time and enhance the localization accuracy. Numerical simulations have revealed that DiSouL is also remarkably robust to calibration errors. The large gain in accuracy comes with higher computational complexity compared to previous existing techniques.

APPENDIX A PROOF OF LEMMA 1

We aim to prove that under A2) and A3), if w > √ L -1, then any estimated location is consistent with L paths (in the sense of Definition 1). Here the point is that it is more costly (in terms of the objective function) to explain an observation as L NLOS angles, than as one position with L associated LOS angles, and this is so exactly when w > √ L -1. Let X and {y} L l= 1 be a solution from [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF] with cost C 1 , and let π 1 be an estimated location (i.e. L l= 1 |x 1l | 2 = 0). Then ẑl can be expressed as ẑl = x 1l a l (θ l (π 1 )) + y 1l a l (ϑ 1l ) + e l , l = 1, . . . , L, [START_REF] Rusek | Scaling up MIMO: Opportunities and challenges with very large arrays[END_REF] where ϑ 1l = θ l (π 1 ) and e l is placeholder of all the terms in the reconstruction ẑl that are not related to the location π 1 or angle ϑ 1l e l = q > 1 x q l a l (θ l (π q )) + m > 1 y m l a l (ϑ m l ) l = 1, . . . , L.

(31) Now, if x 1 0 = L, then π 1 is consistent with L paths in the sense of Definition 1 due to Assumption A3). Here, x 1 = [x 11 • • • x 1L ] T and • 0 is the ℓ 0 -norm which counts the number of estimated elements. Hence, we must prove that w > √ L -1 implies x 1 0 = L. The proof is by contradiction.

Assume that

x 1 0 < L. ( 32
)
This means that the position is consistent with less than L paths. Now we have another competing reconstruction X ′ , {y ′ } L l= 1 and ẑl = (x 1l + y 1l ) .

=y ′ 1 l a l (ϑ 1l ) + e l for l = 1, . . . , L,

with x ′ 1 = 0 and cost C 2 . Since X and {y} L l= 1 are optimal, C 1 ≤ C 2 . Consider now the two assignments [START_REF] Rusek | Scaling up MIMO: Opportunities and challenges with very large arrays[END_REF] and [START_REF] Tropp | Algorithms for simultaneous sparse approximation[END_REF]. Ignoring any common coefficients, the cost of ( 30) is

Define the vector function 1 x whose l-th entry is 1 if x 1l = 0, and 0 otherwise, and denote by x the element-wise absolute value of x 1 , i.e. xl = |x 1l |. Then, x 1 1 = 1 T x x, so from the Cauchy-Schwarz inequality, it follows immediately that

x 1 1 ≤ x 1 0 x 1 2 . ( 37
)
Putting everything together, we find the following contradiction

w x 1 2 (36)
≤ x 1 1 (37)

≤ x 1 0 x 1 2 (32) ≤ √ L -1 x 1 2 (a) < w x 1 2 , (38)
where (a) is due to the fact that w > √ L -1. Hence, w > √ L -1 implies x 1 0 = L.

APPENDIX B PROOF OF LEMMA 2

If no location were found, for each possible location just enough "mass" from each NLOS detected observation could be moved, in a certain way, over to LOS; the new cost cannot exceed the nominal cost if w < √ L. The proof is by contradiction. Assume that w < √ L and that there is no estimated location output by problem [START_REF] Pattipati | A new relaxation algorithm and passive sensor data association[END_REF], so that x q l = 0, ∀q, l. Then, ẑl = m y m l a l (ϑ m l ) l = 1, . . . , L.

Assume without loss of generality that ϑ 1l = θ l (p). By Assumption A3), θ l (p) ∈ Θl , so that y 1l = 0, which leads to the following decomposition ẑl = y 1l a l (ϑ 1l ) + m > 1 y m l a l (ϑ) l = 1, . . . , L.

We consider a competing decomposition X ′ , {y ′ } L l= 1 , for which x ′ q l = 0 for some q, l. In particular, π 1 = p; then ẑl = x ′ 1l a l (θ l (π 1 )) + (y 1lx ′ 1l ) . implying that w √ L -L ≥ 0, which contradicts w < √ L.

  vation, by applying a matched filter (MF) 3 z MF l (t) = T obs 0 s * (τt)z l (τ ) dτ = α l r s (tτ l (p)) a l (θ l (p))

Fig. 1 .

 1 Fig. 1. Example scenario with four base stations. The white circle around each base station excludes the near field region. Each TOA leads to a red disc, the intersection of which is the feasible area. The blue and black lines represent the true bearing lines of the NLOS and LOS paths, respectively.

Algorithm 1 :zl 2 2 > ǫ then 4 :

 124 Direct Localization. 1: set L = L and p = ∅ 2: set ǫ according to (19) 3: if L l= 1 while p = ∅ and L > 1 do 5: set w = L -0.5 6:

  m], [45 m, -45 m], [-45 m, 45 m] and [-45 m, -45 m]. The carrier frequency is 7 GHz and the wavelength 43 mm. Every base station is equipped with a 70-antenna uniform circular array (UCA)

Fig. 2 .

 2 Fig. 2. The left figure plots the LOS and NLOS paths and the right figure plots the bearing lines (BLs). The top left corner base station only receives the LOS path.

Fig. 3 .

 3 Fig. 3. Probability of sub-meter accuracy versus the choice of the weight for the scenario in Fig. 2. Probability estimated by Monte Carlo simulation where the random parameters are the signal strengths and phases. In grey the area comprising the tested values of w satisfying Theorem 1.

Fig. 4 .

 4 Fig. 4. Cumulative density function of the localization error for E/N 0 = 10 dB and B = 30 MHz.

Fig. 5 .

 5 Fig. 5. Probability of sub-meter precision vs. E/N 0 for B = 30 MHz.

Fig. 6 .

 6 Fig.6. Probability of sub-meter precision vs. bandwidth for E/N 0 = 10 dB.

Fig. 7 .

 7 Fig.7. Probability of sub-meter precision vs. number of antennas at each base station for E/N 0 = 10 dB and B = 30 MHz. The diameter of all BSs UCAs grows linearly with the number of antennas according to the formula 0.5π -1 λS l . For 10 antennas the diameter is 7 cm and for 120 antennas is 82 cm.

Fig. 8 .

 8 Fig. 8. Probability of sub-meter precision vs. ray mean arrival time for E/N 0 = 10 dB and B = 30 MHz.

Fig. 9 .

 9 Fig. 9. Contour plot of DiSouL's probability of sub-meter precision vs. gain and phase mismatches on the antennas of all BSs. SNR and bandwidth set to E/N 0 = 10 dB and B = 30 MHz, respectively.

C 1

 1

1 |y

 1 y ′ m l a l (ϑ m l ) l = 1, . . . , L,(41)where y ′ m l = y m l for all l and m > 1. Ignoring common terms, we can associate a cost C 1 and C 2 with (40) and (41), x ′ 1l such that |x ′ 1l | = min l |y 1l | and ∠x ′ 1l = ∠y 1l , and utilize the fact that C 1 ≤ C 2 , we have 1l | -L min l |y 1l |, (44)

Briefly introduced in the conference paper[START_REF] Garcia | Direct localization of emitters using widely spaced sensors in multipath environments[END_REF] but without mathematical justification and for the case of TOA-based localization instead.

Performing the matched filtering requires perfect knowledge of the pulse shape s(t). In practice, if the antennas and hardware have an entirely all-pass (frequency-flat) frequency response, or the pulse shape is slightly changed due to obstacles in the path of the LOS, then the signal-to-noise-ratio will decrease but the number of multipath components will remain the same.

A peak is a local maximum of z NC l (t).

Because on a plane two non-parallel straight lines always intersect, three is the minimum number of BSs needed in order to discriminate LOS from NLOS paths.

While research on massive MIMO has

This work was supported in part by the European Research Council under Grant 258418 (COOPNET) and in part by the EU HIGHTS project (High precision positioning for cooperative ITS applications) under Grant MG-3.5a-2014-636537.