Real-time Range-Angle Estimation and Tag Localization for Multi-static Backscatter Systems

1 Real-time Range-Angle Estimation and T ag Lo calization for Multi-static Bac kscatter Systems T ara Esmaeilb eig ∗ , Kartik P atel ∗ , T raian E. Abrudan, John Kimionis, Eleftherios Kampianakis, Michael Eggleston Abstract—Multi-static bac kscatter netw orks (BNs) are strong candidates for join t communication and localization in the ambien t IoT paradigm for 6G. Enabling real-time lo caliza- tion in large-scale m ulti-static deploymen ts with thousands of devices require highly ecient algorithms for estimating key parameters such as range and angle of arriv al (AoA), and for fusing these parameters into lo cation estimates. W e prop ose t w o low-complexit y algorithms, Join t Range-Angle Clustering (JRA C) and Stage-wise Range-Angle Estimation (SRAE). Both deliver range and angle estimation accuracy compara- ble to FFT- and subspace-based baselines while signicantly reducing the computation. W e then in tro duce tw o real-time lo- calization algorithms that fuse the estimated ranges and A oAs: a maximum-lik eliho o d (ML) metho d solved via gradient searc h and an iterative re-w eighted least squares (IRLS) method. Both ac hiev e lo calization accuracy comparable to ML-based brute force search alb eit with far lo w er complexity . Exp eriments on a real-w orld large-scale m ulti-static testbed with 4 illuminators, 1 multi-an tenna receiver, and 100 tags show that JRAC and SRAE reduce run time b y up to 40 × and IRLS achiev es up to 500 × reduction ov er ML-based brute force search without degrading localization accuracy . The proposed methods ac hiev e 3 m median lo calization error across all 100 tags in a sub-6GHz band with 40 MHz bandwidth. These results demonstrate that m ulti-static range-angle estimation and lo calization algorithms can make real-time, scalable backscatter lo calization practical for next-generation ambien t IoT netw orks. Index T erms—Bac kscatter, integrated comm unications and sensing, Joint range and angle estimation, IoT. I. Introduction The am bien t IoT paradigm envisioned for 6G demands large-scale connectivity and sensing solutions p o w ered by ultra-lo w energy devices. Backscatter comm unication is a strong candidate to realize this vision [ 1 ]–[ 3 ]. While early commercial bac kscatter systems, based on monostatic arc hitectures [ 4 ], were constrained b y short range and lim- ited scalability , the recent adv ances in ambien t and m ulti- static bac kscatter architectures ov ercome these limitations and supp ort wide-area deplo ymen ts while retaining the lo w-p o w er b enets of backscatter devices [ 5 ]–[ 9 ]. Prior art has prov en the feasibility of multi-static and am bien t bac kscatter net w orks (BNs) for communication with fo cus on data transmission and net w ork scalabil- it y [10]–[20] , with few examples exploring localization capabilities in sub-6GHz systems [21]–[23]. Multi-static ∗ T wo authors hav e equal contributions. All au- thors are with Nokia Bell Labs. Email addresses: tara.esmaeilbeig@nokia.com,{kartik.patel, traian.abrudan, ioannis.kimionis, lefteris.kampianakis, michael.eggleston} @nokia- bell-labs.com H 0 H 2 H 1 TX 1 RX T ag F shift TX 2 TX 3 TX 4 Fig. 1. An illustration of a m ulti-static BN comprising 4 TXs, 1 multi-an tenna RX, and a frequency-shifting tag: The tag shifts the carrier signal to an adjacent band to allow the RX process both signals simultaneously , and enable the bistatic range-angle estimation without external synchronization. Com bining range-angle estimates from multiple TX-RX-pairs enable tag localization. lo calization is critical for applications like asset tracking and context-a ware IoT, manifesting a demand for scalable arc hitectures and algorithms that jointly supp ort commu- nication and lo calization. In this work, we introduce a multi-static BN designed to natively supp ort localization of tags. As illustrated in Fig. 1 , the net work consists of multiple single-antenna transmitters (TXs), multi-an tenna receivers (RXs), and frequency-shifting backscatter tags [10]. By shifting and reecting the incident carrier to an adjacent band, each tag enables the RX to pro cess the direct and bac kscattered sig- nals sim ultaneously . This dual-band pro cessing (a) isolates the backscatter signal from strong direct-path interference, and (b) allows time-dierence of arriv al (TDoA)-based bistatic ranging using the direct TX signal as a timing without external sync hronization betw een TXs and RXs. F urthermore, multi-an tenna RXs supp ort angle-of-arriv al (A oA) estimation from the backscatter signal. W e rst develop bistatic range-angle estimation al- gorithms tailored for this arc hitecture. Although range- angle estimation using wideband multi-an tenna channels has been extensively studied [24], [25], scalable solutions suitable for netw orks with h undreds of tags remain an op en c hallenge. T o address this, we propose tw o metho ds: (a) Joint Range-Angle Clustering (JRAC) which is a joint estimator that clusters dominant multipath comp onents in a truncated range-angle heatmap. (b) Stage-wise Range- Angle Estimation (SRAE) which is a decoupled approach that estimates range and angle sequen tially . Compared to con v en tional approaches like 2D FFT [24] and 2D MU- SIC [25] that incur prohibitive computational complexity , JRA C and SRAE achiev e near real-time op erations with- out sacricing accuracy . Bey ond estimation, w e develop tw o real-time multi- static lo calization algorithms that fuse range and angle 2 estimates: (a) an ML-based lo calization solv ed using gradien t ascen t-based algorithm with line search, and (b) an iterative re-weigh ted least squares (IRLS) metho d that deliv ers comparable accuracy with dramatically lo wer computational costs. Finally , w e v alidate our metho ds on a large-scale real- w orld testb ed comprising 4 TXs, 1 RX, and 100 semi- passiv e tags deploy ed in a 20 × 7 m indoor area. De- spite hea vy m ultipath conditions, our prop osed algorithms ac hiev e p ositioning accuracy comparable to baselines while reducing runtime by orders of magnitude. These results demonstrate the practicality of multi-static BNs for inte- grated communication and lo calization in next-generation am bien t IoT netw orks. A. Prior Art 1) Real-world Multi-static Localization Systems: Re- cen t BNs commonly adopt multi-static or am bient bac kscatter architectures for comm unication [12]–[20]. Ho w ev er, only a few works demonstrate real-world lo- calization using m ulti-static architectures in sub-6 GHz bands [21]–[23]. W ork in [21] estimates direction-of- arriv al and lo calizes RFID tags with a multistatic setup that requires a wired shared clo ck betw een radios. The latter, though eective for synchronization of co-located no des, adds an element of deploymen t complexity for sparse, wide-area installations. Additionally , interference from m ultiple transmitters may b ecome signicant when TX and RX antennas are not co-linear, due to bac kscat- tering o ccuring in same channel. In contrast, our sys- tem leverages frequency-shifted backscatter, eliminating the need for external synchronization and reducing TX- induced interference. Alternatively , [22] estimates bistatic ranges using received signal strength (RSS) measurements com bined with particle-lter-based localization, whic h is appealing as a narrowband scheme. Our approach uses wideband signaling for illumination and phase-based ranging, which is more fa vorable in m ultipath environ- men ts [26], due to the full channel sounding charac- teristics. Finally , [23] introduces digital TX-in terference cancellation for in-band tag reections. In our approac h, b y using frequency-shifting tags, we eliminate the need for suc h in terference cancellation sc hemes, maximizing the receiv er dynamic range indep endently for the illumination and bac kscatter signals, resp ectively . 2) Joint Range-Angle Estimation: There are tw o main approac hes for joint range-angle estimation. The rst con- structs a 2D range–angle heatmap and exploits structural assumptions to recov er the Line of Sight (LoS) path. F or example, compressive-sensing metho ds assume sparsity in the range-angle domain to search for the LoS comp onent ecien tly [27], [28]. How ev er, op erating in the sub‑6 GHz bands with less than 40 MHz bandwidth limits the sparsit y in our setup, making such approac hes less eective. The second class of metho ds, such as 2D FFT [24] and 2D MUSIC [25], [29], compute the range–angle heatmap on a quan tized grid and iden tify LoS p eaks through max‑peak detection or iterativ e renemen t, as in [23]. Although accurate, these metho ds incur high computational cost due to dense grid ev aluations. Our prop osed JRAC and SRAE algorithms optimize this class of metho ds by trun- cating the grid searc h and using stage-wise range–angle estimation, resulting in substantial runtime reduction with minimal loss in accuracy . 3) Lo calization T ec hniques: W e mo del errors in AoA measuremen ts using directional statistics [30]. Since the angles are inherently p erio dic, they lie on a circular or spherical manifold rather than in Euclidean space. As a result, widely used statistical distributions are not represen tativ e of the angular errors. T o the best of our kno wledge, directional statistics were rst applied to AoA- based lo calization using v on Mises–Fisher distribution in [31], and using v on Mises distribution in [32], [33]. Sub- sequen tly , maximum likelihoo d (ML) based lo calization with Time-of-Arriv al (T oA) measuremen ts was in tro duced in [34], [35] for active-radio lo calization systems, whic h w as later extended to incorp orate joint range and angle measuremen ts in [36]. In this pap er, we extend this ML- based form ulation to the m ultistatic setting and develop a gradient ascent algorithm with line search for ecien tly estimating the tag p osition. Our second metho d, IRLS, is inspired by the pseudo- linear formulation for bi-static localization introduced in [37]. F ollow-up w ork applied similar ideas to active-radio settings with range-only measuremen ts [38]–[40] and A oA- only measuremen ts [41]. In this pap er, w e extend IRLS to use both range and A oA measuremen ts in a multi- static BN. Drawing from our exp erimental observ ations, w e further rene the IRLS pro cedure to ensure robustness in practical deploymen ts. Lastly , [42], [43] derive Cramer-Rao Low er Bound (CRLB) for m ulti-static lo calization but do not pro vide algorithms to achiev e those limits. Moreov er, temp oral ltering techniques such as Kalman or particle lters can further improv e lo calization performance [22], [44]. In this pap er, we fo cus on ra w lo calization p erformance, without exploring ltering, and assume that the lo calization accu- racy can b e further improv ed using subsequent ltering. B. Summary of Contributions In summary , the contributions of this pap er are as follo ws: • W e introduce a multi-static BN that natively supp orts tag lo calization along with communication, without an y synchronization b etw een TXs and RXs. The system uses the direct carrier signal as a timing reference for bistatic ranging and uses multi-an tenna RXs for AoA estimation. • W e present t w o ecient algorithms, JRAC and SRAE, for joint range–angle estimation in a bistatic bac kscatter system. W e show these methods are computationally ecient without sacricing the es- timation accuracy . • W e develop tw o lo calization algorithms that f use bistatic range and angle. First, we formulate an 3 ML estimator assuming Gaussian range errors and v on Mises–Fisher distributed angular errors, and solv e it via gradient ascen t with line search. Second, w e in tro duce an IRLS-based metho d using pseudo-linear form ulation. W e show that both of these methods signican tly reduce runtime while preserving ov erall lo calization p erformance. • W e build a large-scale exp erimental testb ed and collect tw o datasets: (a) with a 4-tag setup repre- sen ting a simplied scenario, and (b) with a 100- tag setup represen ting a dense practical deploymen t. Using these datasets, we ev aluate our estimation and lo calization metho ds, analyze their computation complexit y (and runtime), and highlight k ey empirical observ ations on the challenging wireless environmen t data. I I. System Mo del A. System Setup The m ulti-static BN consists of N TX TXs and N RX RXs lo cated in D -dimensional space, D ∈ { 2 , 3 } . W e denote the p osition of the i -th TX by p TX ,i , j -th RX by p RX ,j , and a bac kscatter tag by p . W e denote by d ( i,j ) 0 the distance b et w een i -th TX and j - th RX. W e further dene the bistatic range as the distance tra v eled by a signal from i -th TX to j -th RX via the tag. Denoting the bistatic range by d ( i,j ) B , w e hav e d ( i,j ) B = || p − p TX ,i || + || p − p RX ,j || . (1) Let the angular p osition of the TX and the tag, i.e., AoA of the carrier and the backscatter signals, w.r.t. the j - th receiver denoted by ( θ ( i,j ) 0 , ϕ ( i,j ) 0 ) , ( θ ( j ) 2 , ϕ ( j ) 2 ) , where θ ( i,j ) 0 , θ ( j ) 2 ∈ [ − π , π ] denote the azimuthal angles and ϕ ( i,j ) 0 , ϕ ( j ) 2 ∈ [0 , π ] denote the elev ation angles. Then, p − p RX ,j || p − p RX ,j || = h sin ϕ ( j ) 2 cos θ ( j ) 2 , sin ϕ ( j ) 2 sin θ ( j ) 2 , cos ϕ ( j ) 2 i ⊤ . (2) In this work, we are interested in (a) estimating  d ( i,j ) B , θ ( j ) 2 , ϕ ( j ) 2  from the carrier and the bac kscatter signals, and (b) calculating the tag position p using estimated  d ( i,j ) B , θ ( j ) 2 , ϕ ( j ) 2  from all TXs and RXs. T o simplify the notations, w e a v oid sup erscript i, j when the indices of the TX and RX are clear from the context. B. Signal Mo del In the following, w e describ e a bistatic signal mo del for an arbitrary i -th TX and j -th RX. Note that, collectively , the set of all TXs and RXs constitutes a multistatic system, where each individual pair of TX, RX is a bistatic system. W e consider the receiv ers are equipp ed uniformly-spaced N a -elemen t antenna array . Denoting the p osition of i -th elemen t of the antenna array by [ x i , y i , z i ] ⊤ , we dene the ( N a × 1) arra y resp onse vector of the RX as a ( θ , ϕ ) = h exp  ȷ 2 π λ (( x i − x 0 ) cos θ sin ϕ + ( y i − y 0 ) sin θ sin ϕ + ( z i − z 0 ) cos ϕ ) i N a − 1 i =0 . (3) W e consider the TXs transmit an OFDM-based wide- band illumination signal, also called the carrier signal, with the center frequency F c and N s equally-spaced sub carriers spanning the bandwidth B . The tag reects the illumination signal by shifting the center frequency to F c + F shift , and mo dulating a pack et on it. W e call this reected signal the bac kscatter signal. In the following, w e mo del the channels exp erienced by the carrier and the bac kscatter channels. 1) Carrier Channel: W e consider L 0 propagation paths b et w een the TX and RX, and denote the complex- v alued co ecient, the bandwidth-normalized dela y 1 and azim uth/elev ation angles of ℓ -th propagatikon path by α 0 ,ℓ , τ 0 ,ℓ , θ 0 ,ℓ , ϕ 0 ,ℓ , resp ectiv ely . W e extend the single- an tenna c hannel mo del from [10] to write the c hannel frequency resp onse (CFR) of the TX-RX channel, called carrier c hannel, for multi-an tenna RX as H 0 ,n = L 0 − 1 X ℓ =0 α 0 ,ℓ exp  − ȷ 2 π  F c B + n N s  τ 0 ,ℓ  a ⊤ ( θ 0 ,ℓ , ϕ 0 ,ℓ ) , (4) where n denotes the sub carrier index. After down-con version to baseband and concate- nating H 0 ,n for all sub carriers, we hav e H 0 = [ H ⊤ 0 , − N s / 2 , . . . , H ⊤ 0 ,N s / 2 − 1 ] ⊤ ∈ C N s × N a . Then, H 0 satises the factorization H 0 = F ( τ 0 ) D 0 A ⊤ ( θ 0 , ϕ 0 ) , (5) where F ( τ ) = [ f ( τ 0 , 0 , . . . , f ( τ 0 ,L − 1 )] ∈ C N s × L , (6) f ( τ ) = h e ȷ 2 π ( N s / 2) N s τ , . . . , e ȷ − 2 π ( N s / 2 − 1) N s τ i ∈ C N s × 1 , (7) D 0 = Diag ( α 0 , 0 , . . . , α 0 ,L − 1 ) ∈ C L × L , (8) A ( θ 0 , ϕ 0 ) = [ a ( θ 0 , 0 , ϕ 0 , 0 ) , . . . , a ( θ 0 ,L − 1 , ϕ 0 ,L − 1 )] ∈ C N a × L . (9) 2) Backscatter Channel: Similarly , w e consider L 1 propagation paths b et w een the TX and the tag, and denote the complex-v alued co ecient, the bandwidth-normalized delay and AoA of ℓ -th path b y α 1 ,ℓ , τ 1 ,ℓ , θ 1 ,ℓ , ϕ 1 ,ℓ , resp ectiv ely . Then, the TX-T ag CFR is H 1 ,n = L 1 − 1 X ℓ =0 α 1 ,ℓ exp  − ȷ 2 π  F c B + n N s  τ 1 ,ℓ  , (10) where n is the sub carrier index. 1 If the delay is t , then the normalized delay is τ = B t assuming the sampling rate of B . In practice, the RX sampling rate is higher than the bandwidth B . F or the range and angle estimation, how ev er, the sampling rate any more than the bandwidth do es not provide further b enet. Hence, we assume the sampling rate of B . 4 Similarly , we consider L 2 propagation paths b etw een the tag and the RX, and denote the complex-v alued co ecien t, the bandwidth-normalized delay and A oA of ℓ -th path by α 2 ,ℓ , τ 2 ,ℓ , θ 2 ,ℓ , ϕ 2 ,ℓ , resp ectively . F urthermore, the tag shifts the illumination signal by F shift . Therefore, the T ag-RX CFR can b e dened as H 2 ,n = (11) L 2 − 1 X ℓ =0 α 2 ,ℓ exp  − ȷ 2 π  F c + F shift B + n N s  τ 2 ,ℓ  a ⊤ ( θ 2 ,ℓ , ϕ 2 ,ℓ ) . After down conv ersion to baseband, the CFR of the TX-T ag-RX channel, also called the backscatter channel, is dened as H B ,n = H 1 ,n H 2 ,n (12) = L 1 − 1 X ℓ =0 L 2 − 1 X m =0 " α 1 ,ℓ α 2 ,m exp  − ȷ 2 π n N s ( τ 1 ,ℓ + τ 2 ,m )  exp  − ȷ 2 π  F shift B  τ 2 ,m  a ⊤ ( θ 2 ,m , ϕ 2 ,m ) # . (13) Let a bistatic path, indexed b y ℓ ′ = ℓ + L 1 m , denote the total signal path ov er ℓ -th path in the TX-T ag c hannel and m -th path in the T ag-RX channel. Then, the eective channel co ecient, the normalized delay , and the AoA of ℓ ′ -th bistatic path can b e given by β ℓ ′ = α 1 ,ℓ α 2 ,m e − ȷ 2 π F shift τ 2 ,m /B , τ B ,ℓ ′ = τ 1 ,ℓ + τ 2 ,m , and θ B ,ℓ ′ = θ 2 ,m , ϕ B ,ℓ ′ = ϕ 2 ,m . As a result, (13) can b e simplied to H B ,n = L 1 L 2 − 1 X ℓ =0 β ℓ exp  − ȷ 2 π n N s τ B ,ℓ  a ⊤ ( θ B ,ℓ , ϕ B ,ℓ ) . (14) Concatenating H B ,n for all sub carriers, we hav e H B = [ H ⊤ B , − N s / 2 , . . . , H ⊤ B ,N s / 2 − 1 ] ⊤ ∈ C N s × N a . W e dene L = L 1 L 2 as the total n umber of bistatic paths in the bac kscat- ter c hannel. Then, H B satises the factorization H B = F ( τ B ) D B A ⊤ ( θ B , ϕ B ) , (15) where, D B = Diag ( β 0 , . . . , β L − 1 ) ∈ C L × L . In practice, the RX do es not hav e access to the channel information directly , but pro cesses the received signal to get the channel estimates [10]. W e denote an estimate of the c hannel H b y p H . F urthermore, the RX processes N sym OFDM symbols to get a c hannel estimate. Let p H [ k ] denote the estimate of H using k -th OFDM symbol. Therefore, the RX has access to p H 0 [ k ] , p H B [ k ] , ∀ k ∈ { 1 , ..., N sym } to estimate τ B , 0 , θ B , 0 and ϕ B , 0 . W e assume the channels remain static for the duration of N sym sym b ols. In next section, w e discuss the joint range-angle esti- mation algorithms using these c hannel estimates. Without loss of generality , w e subsequently assume the rst path in all c hannels are the LoS path of the c hannels. Therefore, τ i = τ i, 0 , θ i = θ i, 0 , ϕ i = ϕ i, 0 , i ∈ { 0 , 1 , 2 , B } . F urthermore, τ i = B d i /c . I II. Joint Bistatic Range and angle Estimation F or tag lo calization, w e are in terested in estimating the parameters ( τ B , 0 , θ B , 0 , ϕ B , 0 ) of the LoS path in the bac kscatter c hannel corresp onding to an arbitrary TX-RX pair. Notably , the bistatic range (or delay τ B , 0 ) estimation uses time dierence of arriv al (TDoA) metho d with the arriv al of carrier signal as the reference [10]. This k ey design feature enables the ranging without requiring a separate clock sharing infrastructure b etw een TXs and RXs in the net w ork. A ccordingly , let the TDoA estimate b et w een the signal ov er the backscatter channel H B and the signal ov er the carrier channel H 0 b e denoted by ∆ p τ . Then, given the known TX-RX distance d 0 , the estimate of the bistatic range is p d B = c B (∆ p τ + τ 0 , 0 ) . (16) Estimating the TDoA ∆ p τ in volv es estimating the time of arriv al (T oA) o ver the bac kscatter channel H B and the carrier channel H 0 . Each T oA can b e estimated using the channel impulse resp onse (CIR)-based metho d as in v estigated in Section I I I.C of [ 10]. In this work, in addition to estimating the bistatic range, w e also estimate the AoA of the LoS path in the backscat- ter channel. This requires tw o T oA estimates, one from the carrier and one from the bac kscatter channel, and one A oA estimate from the backscatter c hannel. T o achiev e this, we consider tw o approaches: (i) joint estimation of range and angle, and (ii) stage-wise estimation. The join t approac h is employ ed by the baseline algorithms 2D-FFT and 2D-MUSIC describ ed in Section I I I-A and I I I-B, as w ell as the prop osed algorithm JRAC in Section II I-D. The stage-wise approach is used in the prop osed algorithm SRAE in Section I I I-C. F or the purp oses of this work, w e consider a ULA at the RXs and thus, estimate only azimuthal angles. Therefore, in this section, w e assume ϕ B , 0 = π / 2 . Ho w ever, the prop osed metho ds can b e extended to enable estimating ϕ B , 0 with minimal mo dications. Finally , we dene the dela y and angle grid as the ordered sets T =  N s i G τ     i = 0 , . . . , G τ − 1  and (17) Θ =  sin − 1  2 j G θ      j = − G θ / 2 , . . . , G θ / 2 − 1  , (18) resp ectiv ely , suc h that |T | = G τ and | Θ | = G θ . A. 2D FFT Let the 2D sp ectrum of p H i b e S 2DFFT i ( τ , θ ) = N sym X k =1       N s / 2 − 1 X m = − N s / 2 N a − 1 X n =0 [ x H i [ k ]] m + N s / 2+1 ,n +1 e ȷ 2 π ( τ m N s + vn 2 )       , (19) where, τ ∈ T , v = sin θ for θ ∈ Θ , and i ∈ { 0 , B } . W e dene the set of ( τ , θ ) where | S 2DFFT i | has a p eak ab ov e a threshold T 2DFFT min as P 2DFFT i =  ( τ j , θ ℓ ) ∈ T × Θ     | S 2DFFT i ( τ j , θ ℓ ) | > T 2DFFT min , | S 2DFFT i ( τ j , θ ℓ − 1 ) | < | S 2DFFT i ( τ j , θ ℓ ) | > | S 2DFFT i ( τ j , θ ℓ +1 ) | , | S 2DFFT i ( τ j − 1 , θ ℓ ) | < | S 2DFFT i ( τ j , θ ℓ ) | > | S 2DFFT i ( τ j +1 , θ ℓ ) |  . (20) 5 W e then estimate the T oA of the carrier signal p τ 0 , 0 from P 2DFFT 0 as p τ 0 , 0 = arg min τ :( τ ,θ ) ∈P 2DFFT 0 τ . (21) Subsequen tly , we estimate the T oA and the AoA of the bac kscatter signal ( p τ B , 0 , p θ B , 0 ) from P 2DFFT B as p τ B , 0 = min ( τ ,θ ) ∈P 2DFFT B τ > p τ 0 , 0 τ , p θ B , 0 = arg max θ :( p τ B , 0 ,θ ) ∈P 2DFFT B | S 2DFFT B ( p τ B , 0 , θ ) | . (22) W e use the resultant triplet ( p τ 0 , 0 , p τ B , 0 , p θ B , 0 ) to estimate the bistatic range and the AoA as follo ws: p d B = c B ( τ 0 , 0 + p τ B , 0 − p τ 0 , 0 ) , and p θ B = p θ B , 0 . (23) Finally , ( p d B , p θ B ) presents the resultant range-angle pair as- so ciated to the TX-T ag-RX path. The computational com- plexit y of the 2D FFT metho d is O ( G τ G θ log 2 ( G τ G θ )) . B. 2D MUSIC T wo-dimensional MUSIC algorithm is an extension of MUSIC algorithm to jointly estimate T oA and AoA [25], [29]. W e use this algorithm to estimate ( p d B , p θ B ) for a bistatic system. Let the measurement vector p h i [ k ] = v ec ( p H i [ k ]) ∈ C N s N a × 1 b e a stac k of the c hannel estimates of all the subcarriers at all the antennas, for i = { 0 , B } and k = { 1 , . . . , N sym } . Then, from ( 5 ) and (15), w e ha v e p h i [ k ] = vec ( H i [ k ]) + v i [ k ] (24) ( a ) = vec ( F ( τ i ) D i A ⊤ ( θ i )) + v i [ k ] = U ( τ i , θ i ) γ i [ k ] + v i [ k ] , where v i [ k ] denotes c hannel estimation error, ( a ) assumes static channel ov er N sym sym b ols, γ 0 = α 0 , γ B = β B , and U ( τ , θ ) = [ a ( θ 0 ) ⊗ f ( τ 0 ) , . . . , a ( θ L − 1 ) ⊗ f ( τ L − 1 )] . (25) The collected observ ations from all N sym OFDM symbols is giv en by Y i = h p h i [1] , . . . , p h i [ N sym ] i = U ( τ i , θ i )( γ i ⊗ 1 ⊤ N sym ) + V i , (26) where V i = [ v i [1] , . . . , v i [ N sym ]] ∈ C N s N a × N sym . The 2D MUSIC algorithm is based on decomp osition of the sample cov ariance matrix, R Y i = 1 N sym Y i Y H i , into the signal subspace E S and noise subspace E N as R Y i = E S ,i Λ S ,i E H S ,i + E N ,i Λ N ,i E H N ,i , (27) where Λ S ,i and Λ N ,i are diagonal matrices with the eigen v alues of the corresp onding subspaces. Assuming the n um ber of signal space K i , Λ S ,i and E S ,i are determined b y using the K i -largest eigen v alues and asso ciated eigen- v ectors of R Y i , whereas Λ N ,i and E N ,i are determined using the remaining N s N a − K i eigen v alues and the asso ciated eigenv ectors. In practice, the v alue of K i ’s is t ypically unknown and v arious criterion such as Minimum Description Length (MDL), Akaik e Information Criterion (AIC) are used for estimating K i ’s. The orthogonality betw een u ( τ , θ ) = a ( θ ) ⊗ f ( τ ) and the noise subspace is measured as ∥ u ( τ , θ ) H E N ∥ 2 2 , and its in v erse forms the pseudo-sp ectrum of 2D MUSIC as S 2DMUSIC i ( τ , θ ) = ( u ( τ , θ ) H E N ,i E H N ,i u ( τ , θ )) − 1 , (28) for ( τ , θ ) ∈ T × Θ and i ∈ { 0 , B } . Finally , to estimate ( p d B , p θ B ) from S 2DMUSIC i , i ∈ { 0 , B } , the pro cedure described in (20)-(23) is follow ed, albeit with the quantities T 2DMUSIC min and P 2DMUSIC i . The computational complexity of the 2D-MUSIC al- gorithm is dominated by three main comp onents. First, the eigenv alue decomposition (EVD) of the sample co- v ariance matrix R Y i incurs a cost of O ( N 3 a N 3 s ) . Second, the matrix multiplication in (28) inv olving matrices of size N a N s × 1 and N a N s × ( N a N s − K i ) con tributes to the complexity of O  N 2 a N 2 s  . Finally , the pseudo- sp ectrum calculation, which ev aluates the sp ectrum ov er G τ G θ grid p oints, scales the computational complexit y linearly with the num b er of grid points. Therefore, the computational complexity of the 2D MUSIC algorithm is O ( G τ G θ N 2 a N 2 s + N 3 a N 3 s ) . This computational cost may p ose signicant limitations for practical implementations. Therefore, in the following sections, we inv estigate alter- nativ e approaches for joint range and angle estimation that oer improv ed computational eciency . C. Stage-wise Range-Angle Estimation (SRAE) This algorithm op erates in tw o stages: 1) Stage 1 (Dela y estimation): T o estimate the LoS dela y p τ i , we apply the (one-dimensional) MUSIC algorithm in the delay domain. The sample co v ariance matrix is dened as R x H i = 1 N sym N sym X k =1 p H i [ k ] p H H i [ k ] ∈ C N s × N s . (29) The cov ariance matrix is then decomp osed into the signal subspace E S and noise subspace E N using the MDL criterion. The MUSIC sp ectrum is then calculated as S SRAE i ( τ ) = ( f ( τ ) H E N ,i E H N ,i f ( τ )) − 1 . (30) Then, the p eak search is performed on the delay-spectrums of the carrier and backscatter channels, such that p τ 0 , 0 = min  τ j ∈ T     S SRAE 0 ( τ j ) > T SRAE min , S SRAE 0 ( τ j − 1 ) < S SRAE 0 ( τ j ) > S SRAE 0 ( τ j +1 )  , (31) p τ B , 0 = min  τ j ∈ T     τ j > p τ 0 , 0 , S SRAE B ( τ j ) > T SRAE min , S SRAE B ( τ j − 1 ) < S SRAE B ( τ j ) > S SRAE B ( τ j +1 )  . (32) Using the delay estimates p τ 0 , 0 , p τ B , 0 , the bistatic range estimate can b e determined using (23). 2) Stage 2 (AoA estimation): The complex gains of the c hannel taps associated to the LoS path for all an tennas can be retrieved using the channel estimate p H B . These taps associated to the delay p τ B , 0 is giv en by p X B [ k ] = f ( p τ B , 0 ) † p H B [ k ] , k ∈ { 1 , . . . , N sym } , (33) The AoA estimation of the LoS path of the backscatter c hannel, p θ B , is obtained based on the relative phases of p X B , again by applying either b eamforming or subspace-based metho ds. In this instance, we run another round of MUSIC in a spatial domain. The relev ant sample cov ariance R X B = 1 N sym P N sym k =1 X H B [ k ] X B [ k ] , is decomp osed in to the 6 signal subspace E S , B and noise subspace denoted by E N , B . Note that, given there is only single-dela y (and hence, single path) in (33), the size of the signal subspace is 1. The orthogonalit y b et w een a ( θ ) and the noise subspace is measured as ∥ a ( θ ) H E N ∥ 2 2 , the in verse of which forms the MUSIC sp ectrum in spatial domain. This sp ectrum is searc hed for the p eaks as p θ B = argmax θ ∈ Θ 1 a ( θ ) H E N , B E H N , B a ( θ ) . (34) SRAE metho d is equiv alen t to applying t w o one- dimensional MUSIC algorithms. The computational com- plexit y of range estimation is O ( G τ N 3 s ) , while the com- plexit y of the A oA estimation is O ( G θ N 3 a ) . Therefore, the o v erall computational complexity is O ( G τ N 3 s + G θ N 3 a ) . With the SRAE algorithm, the accuracy of delay es- timation aects the accuracy of the estimated taps, and consequen tly , of the A oA estimation. Hence, this stage- wise approac h remains sub-optimal, in fact, irresp ectiv e of the order of estimating the parameters. Thus, w e next in tro duce a joint estimation metho d based on clustering on the joint range-angle heatmap. D. Joint Range-Angle Clustering (JRAC) The JRA C method follows a four-step process: First, similar to SRAE, it iden ties the appropriate search space in delay domain. Second, it estimates a 2D heatmap for the estimated delay space and the entire spatial domain. Third, it clusters v arious channel taps. Finally , it selects the cluster with the low est delay to retrieve the delay and angle corresp onding to the LoS path. In the following, w e formally describ e each of these steps in detail. The outline of our algorithm is given in Algorithm 1 . 1) Step 1 (T runcate the search space in delay-domain): W e rst truncate the p ossible dela y space using 1D MUSIC in the dela y domain. Sp ecically , we use S SRAE i , i ∈ { 0 , B } from (30), and the threshold T JRAC min to determine the truncated dela y space as T i =  τ ∈ T | S SRAE i ( τ ) ≥ T JRAC min  , i ∈ { 0 , B } . (35) 2) Step 2 (The range-angle heatmaps and clipping): W e dene the 2D heatmap asso ciated to the estimated c hannels p H i as S JRAC i ( τ , θ ) = N sym X k =1    f ( τ ) H x H i [ k ] a ∗ ( θ )    2 , τ ∈ T i , θ ∈ Θ , i ∈ { 0 , B } . (36) Note that, this heatmap is similar to the o v ersampled 2D FFT spectrum describ ed in (19), alb eit with a truncated dela y space with reduced computational cost. Hence, it is possible to estimate the bistatic range and angle using equations (20)-(23). How ever, in the following, we in tro duce a clustering step whic h can provide a robust estimate in an heavily multipath environmen t. Let T JRAC max denote a low er clipping threshold and s max ,i = max τ ∈T i ,θ ∈ Θ S JRAC i ( τ , θ ) denote the maximum v alue of the heatmaps for i ∈ { 0 , B } . W e dene a clipped normalized spectrum as ¯ S JRAC i ( τ , θ ) = S JRAC i ( τ , θ ) 1 { S JRAC i ( τ , θ ) ≥ T JRAC max s max ,i } . (37) Notably , ¯ S JRAC i , i ∈ { 0 , B } ha v e islands of clusters. Algorithm 1 Joint Range-Angle Clustering (JRAC) Input p H i [ k ] , τ 0 , 0 , T , Θ , Thresholds T JRAC min and T JRAC max . 1: for i ∈ { 0 , B } do Step 1: T runcate the search space in delay-domain 2: Calculate S SRAE 0 ( τ ) and S SRAE B ( τ ) according to (30). 3: T runcate the delay search space T to get T 0 , T B as T 0 =  τ ∈ T   S SRAE 0 ( τ ) ≥ T min  , T B =  τ ∈ T   S SRAE B ( τ ) ≥ T min  . Step 2: The range-angle heatmaps and clipping 4: Get the 2D heatmap S JRAC i ( τ , θ ) , ∀ τ ∈ T i , θ ∈ Θ using (36). 5: Find the maxima of the heatmap s max ,i = max τ , θ S JRAC i ( τ , θ ) , then clip the S JRAC i such that ¯ S JRAC i ( τ , θ ) = S JRAC i ( τ , θ ) 1 { S JRAC i ( τ , θ ) ≥ T JRAC max s max ,i } . Step 3.a: Clustering in delay-domain 6: Calculate ¯ s i ( τ ) = P | Θ | ℓ =1 ¯ S JRAC i ( τ , θ ℓ ) , ∀ τ ∈ T i , . 7: Determine boundaries of the delay-domain clus- ters by t i,j =      0 , j = 0 or j = |T i | + 1 , 1 , ¯ s i ( τ j )  = 0 , 0 , ¯ s i ( τ j ) = 0 , , j ∈ { 0 , . . . , |T i | +1 } , 8: b i,j = t i,j +1 − t i,j , j ∈ { 0 , . . . , |T i |} . 9: Determine the range-clusters D i,m using (40). Step 3.b: Clustering in angle-domain 10: Calculate ¯ s i ( θ ℓ ) = P |T i | j =1 ¯ s i ( τ j , θ ℓ ) , ∀ θ j ∈ Θ . 11: Repeat Steps 7 - 9 to detect the angle-clusters. Denote n -th cluster as A i,n . Step 4: Estimation 12: Let k -th cluster C i,k = D i,m × A i,n . 13: Determine the local maxima within each cluster as ( p τ i,k , p θ i,k ) = arg max ( j,ℓ ) ∈C i,k ¯ S JRAC i ( τ j , θ ℓ ) . 14: Find the LoS cluster indexed by k ∗ i = arg min k p τ i,k . 15: end for 16: Calculate p τ 0 , 0 = τ 0 ,k ∗ 0 , p τ B , 0 = τ B ,k ∗ B , p θ B , 0 = θ B ,k ∗ B . 17: Rep eat steps 4-9 on S JRAC 0 ( θ, τ ) to estimate the TX-RX delay p τ 0 , 0 . 18: Calculate the LoS bistatic range and angle as p d B = c B ( τ 0 , 0 + p τ B , 0 − p τ 0 , 0 ) and p θ B = p θ B , 0 . Output The estimated bistatic range and angle as ( p d B , p θ B ) . 3) Step 3 (Clustering): W e now identify the b oundaries of the clusters in the delay and angle domain. a) Delay Domain: F or i ∈ { 0 , B } , let ¯ s i ( τ ) denote the the angle-wise aggregated ¯ S JRAC i ( τ , θ ) , such that ¯ s i ( τ ) = | Θ | X ℓ =1 ¯ S JRAC i ( τ , θ ℓ ) , ∀ τ ∈ T i , i ∈ { 0 , B } . (38) W e then determine the non-zero p oin ts t i,j of the range- clusters and the b oundary indicators b i,j in the dela y domain as t i,j =      0 , j = 0 or j = |T i | + 1 , 1 , ¯ s i ( τ j )  = 0 , 0 , ¯ s i ( τ j ) = 0 , , j ∈ { 0 , . . . , |T i | + 1 } , b i,j = t i,j +1 − t i,j , j ∈ { 0 , . . . , |T i |} . (39) Note that for ev ery b i,α = 1 , there is an associated β > α , suc h that b i,β = − 1 and b i,j ′′ = 0 , ∀ α < j ′′ < β . This 7 -20 -15 -10 -5 0 5 10 X [m] -5 0 5 10 15 20 Y [m] Tag RX TX 1 TX 2 TX 3 TX 4 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 Fig. 2. Elliptical p ositioning scenario with bistatic range measure- ments from four TXs and one RX. The tag’s location is determined as the intersection of the four ellipses, each dened by a TX–RX pair as its fo ci and the hyperplane from the estimated tag AoA. indicates the existence of a start and end b oundaries of all dela y-clusters. Let D i dene the set of the tuples ( α m , β m ) containing the index of the start b oundary α m and the end b oundary β m of m -th delay-cluster. W e dene m -th delay-cluster as D i,m = { j ′′ | α m < j ′′ < β m , ∃ ( α m , β m ) ∈ D i } . (40) Note that, ∪ m D i,m = { j | t i,j = 1 } represen ts all indices j with non-zero v alues at ¯ s i ( τ j ) . b) Angle Domain: Similarly , we denote the delay- wise aggregated ¯ s i ( τ , θ ) by ¯ s i ( θ ) , and determine the n -th angle-cluster A i,n using (38)-(40). 4) Step 4 - Estimation: Let the k -th cluster in the heatmap b e dened using the m -th delay-cluster and n -th angle cluster as C i,k = D i,m × A i,n . The algorithm searc hes for the lo cal maxima within each cluster as ( p τ i,k , p θ i,k ) = arg max ( j,ℓ ) ∈C i,k ¯ S JRAC i ( τ j , θ ℓ ) . (41) Then, the cluster asso ciated to the LoS path corresponds to the one with minimum delay . Therefore, k ∗ i = arg min k p τ i,k , (42) p τ 0 , 0 = τ 0 ,k ∗ 0 , p τ B , 0 = τ B ,k ∗ B , p θ B , 0 = θ B ,k ∗ B . (43) Finally , follo wing (23), p d B = c B ( τ 0 , 0 + p τ B , 0 − p τ 0 , 0 ) , and p θ B = p θ B , 0 , (44) w e get ( p d B , p θ B ) , the resultant bistatic range-angle pair asso ciated to the backscatter path. The computational complexity of the JRA C algorithm is dominated by t wo steps: First, the MUSIC algorithm used in Step 1 for truncating the delay domain has complexity O ( G τ N 3 s ) . Second, the computation of determining the heat map in (36) has complexit y O ( G τ G θ N s N a ) . The Step 3 and 4 are one-dimens ional searc hes with the complexity O ( G τ + G θ ) . Therefore, the ov erall computational com- plexit y of the JRAC algorithm is O ( G τ N 3 s + G τ G θ N s N a ) . IV. Positioning In this section, we propose t w o real-time multi-static p ositioning algorithms. As sho wn in Fig. 2 , the m ulti- static positioning scenario with range-angle measurements reduces to nding the intersection of the ellipses, each dened b y a TX–RX pair as its fo ci, along with the h yp erplane determined by the estimated tag AoA. In the most general scenario, only a subset of receiv ers are equipp ed with an tenna arra ys that enable A oA estimation. Therefore, our subsequen t formulation allows for range and angle measuremen ts from dierent receivers to b e exploited separately and do es not assume that b oth range- angle estimates come from the same receiv er. Let p d n denote the bistatic range estimate of n -th measuremen t from the i n -th TX and j n -th RX. Denoting the noise in the range measurements as η r ,n , we can write p d n = ∥ p − p TX ,i n ∥ + ∥ p − p RX ,j n ∥ + η r ,n . (45) W e assume M r n um ber of range estimates for the p osi- tioning. Let p θ m , p ϕ m denote the m -th azim uthal and elev ation A oA measurements at the j m -th RX. Therefore, the unit v ector from the RX to the tag can b e dened as p v m = ( [cos p θ m , sin p θ m ] ⊤ , D = 2 [sin p ϕ m cos p θ m , sin p ϕ m sin p θ m , cos p ϕ m ] ⊤ , D = 3 . (46) W e assume M ∠ n um ber of angle estimates for the po- sitioning. Note that the AoA p θ B in Section II I, whic h is the estimated angle in the an tenna frame, should be con v erted to angles dened in a global reference frame, as depicted in Fig. 5 . After AoA estimation and b efore p ositioning, a transformation matrix should be applied on p v m to p erform conv ersion from the an tenna frame to global frame. W e denote the transformation matrix by Ω and the unit directional vector in global frame will b e p u m = Ω p v m . A. ML-based Positioning in Multi-static BN using Range and angle Measurements W e rst introduce an ML-based p ositioning algorithm for m ulti-static BNs inspired from [36]. 1) ML F orm ulation: As in (45), we assume that n -th range measurement is corrupted by zero-mean Gaussian noise with v ariance σ n , i.e., η r ,n ∼ N (0 , σ n ) . Therefore, the log lik eliho od function of the tag p osition given the n - th range estimate, after ignoring additive constants w.r.t. p , is given by [34] ℓ ( r ) n ( p ) = −  p d n − ∥ p TX ,i n − p ∥ − ∥ p RX ,j n − p ∥  2 2 σ 2 n . (47) Similarly , we assume that m -th angle measurement p θ m is corrupted suc h that the estimated directional estimate p u m follo ws von Mises-Fisher distribution with the mean direction u m and the concentration parameter κ m . Hence, p u m ∼ VMF ( u m , κ m ) , and (48) P ( p u m | u m , κ m ) = κ m 4 π sinh ( κ m ) exp ( κ m p u ⊤ m u m ) . (49) Therefore, the log likelihoo d function of the tag p osition giv en the angle estimate p θ m , after ignoring additive constan ts w.r.t. p , is given by ℓ ( ∠ ) m ( p ) = log ( P [ p u m | p ]) = κ m p u ⊤ m p − p RX ,j m || p − p RX ,j m || , (50) where w e replace u m b y ( p − p RX ,j m ) / || p − p RX ,j m || . Remark 1. W e adopt a directional statistics approac h b ecause the natural parameter space of angles is not an Euclidean space, but a sphere. Therefore, a distribution 8 on the angles must account for the implicit p erio dicit y of angles and must hav e a natural supp ort ov er the unit sphere. The v on Mises-Fisher distribution is the analog of the circular bi-v ariate normal distribution from the Euclidean space to the unit sphere. In fact, vMF distri- bution is equiv alent to a bi-v ariate Gaussian distribution constrained onto a unit sphere in a D -dimensional space. The mean of such vMF distribution is the unit vector of the mean of the asso ciated Gaussian distribution, and the reliabilit y metric κ m is asso ciated to the v ariance of the Gaussian distribution. Refer to [30], [31] for more details. The log lik eliho o d function L ( p ) of the tag p osition giv en M r -range and M ∠ -angle measuremen ts is given by L ( p ) = M r X n =1 ℓ ( r ) n ( p ) + M ∠ X m =1 ℓ ( ∠ ) m ( p ) , = M r X n =1 −  p d n − ∥ p TX ,i n − p ∥ − ∥ p RX ,j n − p ∥  2 2 σ 2 n + M ∠ X m =1 κ m p u ⊤ m p − p RX ,j m || p − p RX ,j m || . (51) The tag’s p osition can b e estimated by maximizing the joint log likelihoo d function, i.e., p p ML = argmax p ∈ R D L ( p ) . (52) Since the likelihoo d function is not a con v ex function, a brute force approac h using the grid search can b e emplo y ed. This approach includes (1) quantizing the region of interest into grid p oints with resolution ε ML , (2) ev aluating the log lik eliho o d function on eac h p oint in the grid, and (3) selecting the grid p oint that maximizes the the lik elihoo d as a lo cation estimate. While the grid searc h approach is guaranteed to pro vide a global optimal solution (up to the grid resolution), its computational cost increases inv ersely with ε ML as O ( D ( M r + M ∠ ) ε − D ML ) , where O ( D ( M r + M ∠ )) is the cost of ev aluating the log- lik eliho od function. F or a practical deploymen t, w e are interested in lo w complexit y algorithms for real-time op erations. Therefore, w e introduce a gradient ascen t-based heuristic algorithm to ecien tly nd the solution of (52). 2) ML Gradien t Ascen t with Line Search: T o perform gradien t ascent, we rst obtain the gradients of L ( p ) as ∇L ( p ) ≜ d d p L ( p ) (53) = − 1 2 M r X n =1 1 σ 2 n "  p d n − ∥ p TX ,i n − p ∥ − ∥ p RX ,j n − p ∥   p TX ,i n − p ∥ p TX ,i n − p ∥ + p RX ,j n − p ∥ p RX ,j n − p ∥  # (54) + M ∠ X m =1 κ m  I ∥ p − p RX ,j m ∥ − ( p − p RX ,j m )( p − p RX ,j m ) ⊤ ∥ p − p RX ,j m ∥ 3  p u m . (55) Since the log-likelihoo d function is not a con v ex func- tion, the gradient ascent algorithm is not guaranteed to pro vide the optimal solution. Therefore, we introduce tw o Algorithm 2 ML-based Gradient Ascent with Line Search Input { p d n } M r n =1 , { p u m } M ∠ m =1 , S , K 1: Initializations { p p ML ,i (0) } 5 M r i =1 2: for i = 1 to 5 M r do 3: k ← 0 4: repeat 5: µ ∗ k ← arg max µ k ∈S L ( p p ML ,i ( k ) + µ k ∇L ( p p ML ,i ( k ))) 6: p p ML ,i ( k + 1) ← p p ML ,i ( k ) + µ ∗ k ∇L ( p p ML ,i ( k )) 7: k ← k + 1 8: un til ∥ p p ML ,i ( k ) − p p ML ,i ( k − 1) ∥ ≤ ε ML or k = K ML 9: end for 10: i ∗ ← arg max i ∈{ 1 ,..., 5 M r } L ( p p ML ,i ( k )) Output p p ML ← p p ML ,i ∗ ( k ) features in the gradient ascen t to improv e the p ossibility of reac hing to the global maximum. (a) Multiple initializations: The gradien t ascen t iterations are initialized at the cen ter, t w o vertices and t w o co- v ertices of the ellipse asso ciated to each TX-RX pair (see Fig. 2 ). This results in 5 M r initialization p oints denoted b y p p ML ,i (0) , i = 1 , . . . , 5 M r . An indep endent instance of the gradient ascen t algorithm is executed for each initial p oin t. (b) Multiple step-sizes in each gradient ascent iteration: F or each instance of the gradient ascen t, the estimate is up dated at each iterations by selecting a step-size from a xed set S whic h maximizes the log-likelihoo d function. Sp ecically , p p ML ,i ( k + 1) = p p ML ,i ( k ) + µ ∗ k ∇L ( p p ML ,i ( k )) , (56) where, µ ∗ k = arg max µ k ∈S L ( p p ML ,i ( k ) + µ k ∇L ( p p ML ,i ( k ))) . (57) Dynamically selecting the step-size in eac h iteration of gradien t ascent is also called gradien t ascent with line searc h. At the end of eac h instance of the gradient ascent, the p osition that maximizes the likelihoo d is selected as the lo cation estimate. As for the stopping criterion, w e either use ε ML as a predened p osition accuracy threshold or a maxim um num ber of iterations K ML . W e summarize this algorithm, subsequently called Gradien t Ascent with Line Searc h, in Algorithm 2 . Giv en 5 M r initializations, |S | ev aluations of the log- lik eliho od functions in each iteration, and O ( D ( M r + M ∠ )) complexit y for the ev aluation of the log-likelihoo d function, the total computational cost of the prop osed ML-based gradien t ascent with line searc h is given b y O ( D ( M r + M ∠ ) M r |S | K ML ) . Note that, compared to ML-based grid searc h, the gradien t search reduces the complexit y by a v oiding the exp onen tial dependence on the num ber of dimensions D . This reduction in the complexity comes as a cost of the p ossibilit y of a sub-optimal solution and an additional linear dependence on M r . B. Iterative Re-weigh ted Least Squares (IRLS) In this section, we prop ose an iterative p ositioning approac h called IRLS. While sev eral prior w orks consider the IRLS algorithm [37], [39], [41], they primarily fo cus on either range- or angle-only position estimation. W e instead 9 in tro duce the IRLS method that uses b oth range and angle estimates for multi-static p ositioning. A ccordingly , w e rst formulate a pseudo-linear problem for multi-static p ositioning using range and angle measurements. W e then in tro duce iterative w eigh t updates, and nally , summarize the algorithm and discuss its computational complexity . 1) Pseudo-linear F orm ulation [37]: Assuming the noise η r ,n ≪ p d n , b y re-arranging and squaring (45), w e get 1 φ r ,n  2( p RX ,j n − p TX ,i n ) ⊤ 2 p d n   p ∥ p − p RX ,j n ∥  = 1 . (58) where φ r ,n = p d 2 n − ∥ p TX ,i n − p RX ,j n ∥ 2 + 2( p RX ,j n − p TX ,i n ) ⊤ p RX ,j n . (59) Collecting (58) for all M r range measurements, we dene Φ r = [ Φ ⊤ r , 1 , Φ ⊤ r , 2 , . . . , Φ ⊤ r ,M r ] ⊤ ∈ R M r × ( D + N RX ) , (60) Φ r ,n = 1 φ r ,n [2( p RX ,j n − p TX ,i n ) ⊤ , 0 ⊤ j n − 1 , 2 p d n , 0 ⊤ N RX − j n ] . (61) Similarly , assuming lo w angle measuremen t errors, w e ha v e  I D − p u m   p ∥ p − p RX ,j m ∥  = p RX ,j m ≜ φ ∠ ,m . (62) Collecting ( 62) for all M ∠ angle measurements, we dene Φ ∠ = [ Φ ⊤ ∠ , 1 , . . . , Φ ⊤ ∠ ,M ∠ ] ⊤ ∈ R DM ∠ × ( D + N RX ) , (63) Φ ∠ ,m = [ I D , 0 D × ( j m − 1) , − p u m , 0 ⊤ D × N RX − j m ] , (64) φ ∠ = [ φ ⊤ ∠ , 1 , . . . , φ ⊤ ∠ ,M ∠ ] ⊤ ∈ R DM ∠ × 1 . (65) Finally , let Υ =  p ⊤ , ∥ p − p RX , 1 ∥ , . . . , ∥ p − p RX ,N RX ∥  ⊤ , (66) Φ =  Φ ⊤ r , Φ ⊤ ∠  ⊤ ∈ R ( M r + DM ∠ ) × ( D + N RX ) , (67) φ =  1 ⊤ M r , φ ⊤ ∠  ⊤ ∈ R ( M r + DM ∠ ) × 1 . (68) F rom (58) and (62), (66)-(68) can b e collectively written as ΦΥ = φ . (69) T o estimate the tag position, w e need to estimate Υ using a least square solution given Φ and φ . This provides a pseudo-linear 2 form ulation of the m ulti-static p osition- ing problem, given the range and angle measuremen ts. F urthermore, this formulation can also b e extended to include the measuremen ts from multiple metho ds simply b y stacking more rows in the matrix Φ and φ . Remark 2. In addition to using both range and angle measuremen ts, our IRLS form ulation diers from the prior w ork in another w a y . All rows corresp onding to the range measuremen ts in (69) are scaled such that their resp ective elemen ts in φ is 1. F rom (58), note that φ r ,n dep ends on the locations of the TXs and RXs, and p d n . While most prior work fo cused on simulations which resulted in comparable φ r ,n across all measurements, our exp erience with the real-w orld prototype show ed that, in practice, φ r ,n v aries drastically across measuremen ts. Hence, solving (69) without the scaling by φ r ,n w ould give articially high importance to the measuremen ts with high φ r ,n . This problem is further exacerbated if n -th range measurement is an outlier resulting in high p d 2 n , which w ould add higher w eigh t to an outlier measurement! Therefore, w e tac kle this problem (highlighted by our exp erimental setup) b y 2 (69) is a pseudo-linear formulation b ecause the parameters dened in Υ are not indep endent v ariables, but hav e non-linear dependency b et ween them. simply scaling (58) by φ r ,n suc h that the right hand side is one. While a least squares solution solv es (69), from a practical persp ective, there are tw o key challenges: • Errors in the measurements: (69) do es not account for the p ossible errors in the range and angle measurements due to signal noise and m ultipaths in the propagation en vironmen t. • V alidity of estimated Υ : The least square solution of (69) do es not enforce the dep endency within { Υ i } D + N RX i =1 . Sp ecically , the solution p Υ is only v alid if ∥ [ p Υ ] 1: D − p RX ,j ∥ = [ p Υ ] D + j , ∀ j ∈ { 1 , . . . , N RX } . (70) In the following, w e introduce tw o metho ds to tackle these t w o challenges. 2) Iterative Re-W eigh ted Least Square: The formu- lation in (69) assumes all measurements to ha v e same imp ortance. How ever, an error in one measurement ma y distort the estimate. Therefore, we introduce a measure of reliabilit y for each measurement. W e dene a weigh t vector w ∈ R M r + D R M ∠ as the w eigh ts assigned to each row of (69). Then, with W = Diag { w } , we solve p Υ = min Υ ( φ − ΦΥ ) ⊤ W ( φ − ΦΥ ) , (71) whic h results in p Υ =  √ WΦ  † √ W φ . (72) The position estimate is then given b y p p IRLS = [ p Υ ] 1: D . (73) T o estimate the p osition p p IRLS robustly , our goal is to select the w eigh ts w r and w ∠ suc h that the outliers are assigned lo wer weigh t. F or that, w e use an iterative approac h as follows: 1) Initialize the weigh t vector by w = 1 M r + D R M ∠ . 2) Estimate p Υ using (72). 3) Estimate the residual errors e = | ΦΥ − φ | . 4) Dene w suc h that w i = ( e 2 i + ε IRLS ) − 1 . 5) Rep eat 2-4 for K num b er of iterations, where ε IRLS prev en ts excessiv ely large w eights and ensures stabilit y . W e emphasize here that while we set same initial w eigh ts for all measurements, in practice, there are other approac hes using measurement metrics. F or instance, the measuremen ts from the tag pack ets with higher signal strength quality or lo w er bit error rate can be assigned higher initial weigh ts. How ev er, given that IRLS do es not guaran tee an optimal solution, it is not p ossible to denitiv ely conclude the b est c hoice for the initial w eights. 3) Ensuring V alidity of the Estimate: The solution of (72) assumes that the elements of Υ are independent. Ho w ev er, the estimated p Υ is v alid only if (70) is satised. T o ensure that this constraint is satised, we use a second stage least squares solution [37], [39] to pro ject ˆ Υ on to feasible cone of (70). Let η IRLS represen t the error betw een the estimated p Υ and the true Υ i.e. p Υ = Υ + η IRLS . W e dene Υ 2 = p 2 ⊙ p 2 as the unkno wn of the second stage least squares and corresp ondingly p 2 as the tag p osition estimate after pro jection. Then, 10 Υ 2 = [ p Υ ] 1: D ⊙ [ p Υ ] 1: D − 2[ η IRLS ] 1: D ⊙ [ p Υ ] 1: D . (74) Squaring (70) yields, 1 ⊤ Υ 2 − [ p Υ ] 2 D + j − 2 p ⊤ RX ,j [ p Υ ] 1: D + || p RX ,j || 2 = − 2[ p Υ ] D + j [ η IRLS ] D + j − 2 p ⊤ RX ,j [ η IRLS ] 1: D , (75) for j = { 1 , . . . , N RX } . In (74)-(75), we use Υ = p Υ − η IRLS and elliminated the second-order error terms. Conse- quen tly , (74)-(75) are recast in matrix form as Φ 2 Υ 2 = φ 2 , where Φ 2 = [ I D , 1 D , . . . , 1 D ] ⊤ ∈ R ( D + N RX ) × D , (76) φ 2 = p Υ ⊙ p Υ +        0 D 2 p ⊤ RX , 1 [ p Υ ] 1: D − || p RX , 1 || 2 2 p ⊤ RX , 2 [ p Υ ] 1: D − || p RX , 2 || 2 . . . 2 p ⊤ RX ,N RX [ p Υ ] 1: D − || p RX ,N RX || 2 .        (77) Then, w e nd p 2 suc h that p p 2 = arg min p 2 ∈ R D || Φ 2 Υ 2 − φ 2 || 2 , (78) whic h results in the solution p Υ 2 = p p 2 ⊙ p p 2 = Φ † 2 φ 2 . (79) Therefore, the nal position estimate from the IRLS algorithm can b e given as p p IRLS = sign ([ p Υ ] 1: D ) ⊙ q Φ † 2 φ 2 . (80) Note that (80) impro ves up on the estimate (73) b y fullling the constraint introduced in (70). 4) Algorithm Complexity: The computation for IRLS algorithm is dominated b y the pseudo-in v erse op eration dened on line 3 of the algorithm. Given that the dimen- sions of Φ is ( M r + D M ∠ ) × ( D + N RX ) , the complexit y of the pseudo-inv erse operation is O ( D ( M r + D M ∠ ) 2 ) . Therefore, the total complexity of the IRLS algorithm is giv en by O ( K D ( M r + DM ∠ ) 2 ) . V. Numerical Results In this section, we present the simulated p erformance of the prop osed estimation and lo calization algorithms. A. Simulation Setup F our TXs are lo cated at the coordinates (0 , 0) , (0 , 6) , ( − 18 , 6) , and ( − 18 , 0) m within the ro om. The RX, equipp ed with an an tenna array with N a = 4 half- w a v elength–spaced elements, is p ositioned at ( − 9 , 3) m. 25 tags are placed uniformly in the region ( x, y ) ∈ [ − 18 , 0] × [0 , 6] m. All TXs, RX and tags are assumed to hav e same height, making the entire deploymen t a 2D deplo ymen t. Scatterers are uniformly distributed in the same region as the tag p ositions. W e assume L 1 , L 2 = 3 num b er of scatterers in the environmen t. F urther, we generate 50 OFDM symbols, with 40 sub carriers uniformly spaced at 1 MHz, according to (24). W e further contaminate the OFDM symbols by Gaussian noise v ( n ) ∼ N (0 , σ 2 v ) such that the SNR = ∥ H ( n ) B ,n ∥ σ 2 v is 5 dB. B. Range-Angle Estimation Performance W e ev aluate joint range–angle estimation metho ds from Section I I I: 2D FFT, 2D MUSIC, SRAE with T SRAE min = 0 . 5 , and JRA C with T JRAC min = 0 . 2 and T JRAC max = 0 . 6 . W e also add a range-only estimation using IR First [ 10] as a baseline to em ulate a single-antenna RX. F or all algorithms, w e use same sets T and Θ with G τ = 4096 , G θ = 128 . This ensures that p erformance and runtime of prop osed algorithms and IR First are fairly compared. Fig. 3 compares the performance of these ve algorithms in terms of the absolute range estimation error | p d − d | and the absolute angle estimation error | p θ − θ | . Fig. 3a shows that the range-angle estimation algorithms outp erforms the range-only baseline. This indicates the b enet of separating the m ultipath comp onents along the spatial dimensions to improv e the range estimates. Among the range-angle estimation algorithms, both Fig. 3a and 3b sho ws (a) JRAC outperforming 2D FFT and 2D MUSIC, and (b) SRAE under-p erforming 2D FFT and 2D MUSIC. JRA C outp erforms due to the clustering step which considers a multiple paths as a unit instead of ev aluating individual paths as done by 2D FFT and 2D MUSIC. F urthermore, SRAE under-p erforms the joint metho ds due its stage-wise architecture. In theory , since delay and angle are coupled in multipath environmen ts, joint algorithms can separate them more accurately . C. Positioning Performance The b enet of p ositioning based on joint range and angle measuremen ts from synthetically generated dataset, com- pared to range-only measurements, is illustrated in Fig. 3c. In a t wo-dimensional positioning scenario, range-only lo calization requires a minimum of three ellipses corre- sp onding to three TXs to estimate a tag p osition. In con trast, when the RX is equipp ed with an an tenna array and the A oA estimates are av ailable, the tag p osition can b e determined using as few as one range and one angle measuremen ts from a single TX. Fig. 3c further sho ws that ev en with a sucient n um b er of TXs, equipping the RX with a ULA to enable A oA estimation signican tly im- pro v es p ositioning accuracy . Sp ecically , with three TXs, range-only measurements result in an av erage p ositioning error of approximately 6 m, which is reduced to ab out 2 m when A oA estimation at the RX is feasible. VI. Exp erimental Results W e conduct the experimental ev aluation on a custom- built wireless testbed deplo yed in an enclosed lab space (metal walls). The testb ed spans an area of 20 m × 7 m and includes four TXs, a single m ulti-channel RX, and 100 custom-designed semi-passive bac kscatter tags distributed across the area. The TXs and RX op erate without an y clo c k sync hronization, reecting a realistic m ulti-static arc hitecture without clo ck distribution. The exp erimental setup is given in Fig. 4 . Eac h TX uses an ADR V9361-Z7035 chip connected to a dip ole antenna and contin uously broadcasts an OFDM wa v eform with 39.36 MHz bandwidth centered at 11 0 2 4 6 8 10 12 Range Error (m) 0 0.2 0.4 0.6 0.8 1 CDF SRAE JRAC 2D MUSIC 2D FFT IR First (a) CDF of range estimation er- ror 0 5 10 15 20 Absolute AoA Error (deg) 0 0.2 0.4 0.6 0.8 1 CDF SRAE JRAC 2D MUSIC 2D FFT (b) CDF of angle estimation er- ror 0 2 4 6 8 10 Positioning error (m) 0 0.2 0.4 0.6 0.8 1 CDF Range Only, 3 TXs Range Only 4 TXs Joint Range & AoA, 1 TX Joint Range & AoA, 2 TXs Joint Range & AoA, 3 TXs Joint Range & AoA, 4 TXs (c) CDF of localization error Fig. 3. (a)-(b) Range and angle estimation errors using simulation data: Range-AoA estimation metho ds hav e higher ranging accuracy compared to the range-only metho d IR-First. Joint metho ds such as JRAC achiev es higher ranging and angle estimation accuracy compared to the stage-wise metho d, SRAE. (c) Using AoA for lo calization with v arious num b er of TXs: W e compare p ositioning accuracy using range only and joint range-angle measurements from the synthetic data. W e consider SNR = 5 dB, L 1 = L 2 = 3 scatterers, 1D MUSIC for the range estimation and 2D MUSIC for joint range and angle estimation. Localization uses ML-based Gradient Ascent with Line Search. Fig. 4. The exp erimental setup of a large-scale BN consisting of 4 TXs, 1 RX with 4-element antenna arra y , and 100 tags in an enclosed lab. The walls with metal surface and the ground introduce signicant multipath in the channel. -18 -16 -14 -12 -10 -8 -6 -4 -2 0 X (m) -1 0 1 2 3 4 5 6 7 8 Y (m) 1 2 3 4 Illuminators Receivers 100-tag Setup 4-tag Setup Antenna Frame Global Frame Y X Y X Fig. 5. Schematic of the exp erimental setup consisting 4 TXs, 1 RX with 4-element an tenna array , and 100 tags. The positive-X direction in the antenna frame denotes the broadside direction of the receiver antenna array . 897.5 MHz. The OFDM wa veform contains 41 sub carriers spaced at 960 kHz, following the design principles in [10]. The backscatter tags hav e the semi-passive architecture describ ed in [10]. Each tag transmits one 100-bit pack et of duration 4.8 ms at a uniformly random time within ev ery 300 ms window [11]. T ags mo dulate the inciden t OFDM carrier using a square-wa ve switc hing mec hanism and apply a frequency shift of 45 MHz to separate the backscatter signal from the carrier, which preven ts in terference with the direct signal from the TX. The receiver is built around a USRP X440 with eight RF c hannels. A 4-element half-wa v elength-spaced custom ULA feeds into 4 RF diplexers, which split each antenna feed into tw o bands containing: (a) carrier signal from the direct TX-RX path, and (b) backscatter signal from the TX-tag-RX path. This conguration pro duces eigh t input streams to the USRP . Since USRP X440 lacks internal amplication, w e add external LNAs with 20 dB gain eac h, which b oosts the SNR of weak backscatter signals. The USRP samples all eigh t streams at an ADC rate of 61.44 MSps and sends them to the host computer for pro cessing. I/Q pro cessing follows the steps describ ed in [10]. The RX estimates CFRs for all streams, estimates the phase and timing references from the four carrier channels, and computes the range and angle estimates from the four bac kscatter channel estimates of each detected tag. Then, the estimated tag direction in the an tenna frame is transformed to the world frame as explained in (46). F or the top ology illustrated in Fig. 5 , the transformation matrix is Ω =  0 1 1 0  . Finally , the range and angle estimates for eac h tag using the carrier signals from all 4 TXs are used to estimate the p osition of the tags. W e ev aluate t w o scenarios as shown in Fig. 5 : • 4-tag setup: This setup presents an ideal condition with high SNR and few collisions. W e k eep four tags close to the RX on and capture 1.2 seconds of con tin uous samples. • 100-tag setup: This setup emulates a realistic de- plo ymen t with low SNR and dense deplo ymen t. W e k eep all 100 tags on and capture 15 seconds of con tin uous samples to ensure sucient observ ations for ev aluation. In b oth scenarios, the TXs transmit con tinuously and the RX records raw I/Q samples for pro cessing. W e ev aluate four range–angle estimation algorithms describ ed in Section I I I: (a) 2D FFT, (b) 2D MUSIC, (c) SRAE with T SRAE min = 0 . 5 , and (d) JRAC with T JRAC min = 0 . 2 and T JRAC max = 0 . 6 . Similar to Fig. 3a, we also add a range-only estimation using IR First [10] as a baseline to emulate a single-an tenna RX setup. F or all algorithms, we use same sets T and Θ with G τ = 4096 and G θ = 128 . This ensures that p erformance and runtime of proposed algorithms and 12 IR First are fairly compared. In Fig. 6a, w e sho w the CDF of absolute ranging error | p d − d | for the 4-tag scenario. The 2D FFT and IR First algorithms maintain p erformance similar to the sim ulated results in Fig. 3a. In contrast, the 2D MUSIC, SRAE, and JRAC metho ds exhibit noticeable degradation in ranging accuracy . This degradation o ccurs b ecause these algorithms rely on a MUSIC-based step for delay estimation. MUSIC assumes a well-dened signal and noise subspace, which b ecomes unreliable in real-world m ultipath environmen ts. Coherent multipath comp onents and uncertaint y in the num ber of subspaces signicantly reduce the robustness of MUSIC for dela y estimation, as noted in prior studies [45], [46]. In Fig. 6b, w e presen t the CDF of absolute A oA estimation error. Unlike the range results, MUSIC-based algorithms outperform FFT-based metho ds for angle estimation. This impro vemen t occurs b ecause, at high SNR, spatial subspaces are more separable than dela y subspaces, even in multipath en vironmen ts. The ULA pro vides spatial diversit y , which allows MUSIC to exploit subspace decomp osition for high-resolution angle estima- tion. In con trast, FFT-based metho ds suer from coarse grid resolution and sp ectral leakage, which limit angular accuracy . Therefore, while MUSIC struggles with delay estimation, it remains eective for AoA estimation under high-SNR conditions. In Fig. 6d and 6e, w e presen t the CDF of absolute ranging and A oA errors in 100-tag setup. While the relativ e p erformance of v arious metho ds remain similar to that of the 4-tag setup, the ov erall p erformance of the estimation algorithms drops substantially . This is due to three p ossible reasons: (a) low SNR in 100-tag setup, (b) additional multipath introduced by the tags when not transmitting a pack et. W e observe that CIRs of the carrier and backscatter channels in the 100-tag setup are signican tly more cluttered than in the 4-tag setup for the same TX, RX, and the tag. This sho ws that there are more m ultipaths in the en vironmen t with 100 tags, indicating that when the tags are on and not mo dulating a pac k et, they indeed act as (unmo dulated) reectors, due to the semi-passiv e arc hitecture. Unlik e passiv e tags where the tag antenna is connected to a matched load (and therefore absorbing inciden t RF signals) when not mo dulating, semi-passive tags hav e their antenna alwa ys connected to a mismatched load [47]. Therefore, even when the tags are not mo dulating, the signals are reected with a phase shift, resulting in additional multipaths. A p oten tial resolution w ould b e the use of multi-state RF fron t-ends to “mute” the tags that are not transmitting pac k ets, b y introducing a matched or low-reection state (e.g. [48]). F urthermore, in Fig. 6e, we also observ e the longer tail in the angle estimation error. The longer tails denote the fundamental limit of the angle estimation in the end-re region of the ULA, in which many tags are lo cated (refer to Fig. 5 ). T able I presents the av erage runtime required to com- pute a range and angle estimate from each tag pack et T ABLE I Computational complexity and R untime of Joint Range and Angle Estimation Metho ds for N sym = 98 OFDM symbols. Estimation Method Computational Complexity A verage Run time (ms) SRAE O ( G τ N 3 s + G θ N 3 a ) 55 JRAC O ( G τ N 3 s + G τ G θ N s N a ) 475 2D MUSIC O ( G τ G θ N 2 a N 2 s + N 3 a N 3 s ) 10380 2D FFT O ( G τ G θ log( G τ G θ )) 2038 T ABLE I I Run time of Positioning Metho ds Positioning Method Computational Complexity Run time (ms) 4-tag 100-tag ML Grid Search O ( D ( M r + M ∠ ) ε − D ML ) 4670 22802 ML Gradient Ascent with Line Search O ( D ( M r + M ∠ ) M r |S | K ML ) 500 15818 IRLS O ( D K ( M r + DM ∠ ) 2 ) 38 51 using the four algorithms. The SRAE and JRAC metho ds outp erform the others by orders of magnitude, ac hieving near real-time pro cessing capability . In con trast, 2D FFT and 2D MUSIC incur signicantly higher computational costs due to their reliance on complex FFT and exp ensive grid searc hes, respectively . This dierence in runtime b ecomes critical for scalabilit y in dense deplo ymen ts. While the relative accuracy and the runtime of dierent range-angle estimation algorithms v aries, ultimately the metric of interest is the ov erall p ositioning quality using the respective range and angle estimates. W e ev aluate the p ositioning performance of tw o opti- mization methods describ ed in Section IV: (a) ML-based Gradien t Ascent with Line Searc h prop osed in Algorithm 2 with σ n , κ n = 1 and the line searc h step sizes S = { 0 . 001 , 0 . 01 , 0 . 1 , 1 } . (b) IRLS prop osed in Section IV-B with ε IRLS = 10 − 6 . The analysis uses range-angle esti- mates obtained from four algorithms: JRAC, SRAE, 2D MUSIC, and 2D FFT.Fig. 6c and 6f shows the CDF of p ositioning error for the 4-tag and 100-tag datasets. The 4-tag dataset achiev es substan tially better positioning accuracy than the 100-tag dataset. This impro v emen t is exp ected because the 4-tag scenario represents a simplied en vironmen t with high SNR, whereas the 100-tag scenario in tro duces low er SNR and heavy m ultipath, which degrade lo calization accuracy . Among the range-angle estimation algorithms, 2D FFT combined with ML-based positioning ac hiev es the low est p ositioning error. This result aligns with the robustness of 2D FFT on real-world data and the optimalit y of ML-based estimation. Interestingly , the tw o lo w-complexit y algorithms proposed in this pap er–JRAC and SRAE–deliver positioning accuracy comparable to 2D FFT when paired with IRLS. This nding demonstrates that reducing computational complexit y in range-angle estimation do es not necessarily compromise positioning qualit y , making JRAC and SRAE practical for large-scale deplo ymen ts. T able II rep orts the av erage run time for computing p osition estimates using ML Gradient Ascent with Line Searc h and IRLS for b oth datasets. It also includes a brute-force baseline that searches the ML ob jective o v er a grid with 5 cm granularit y . The results sho w that ML- based metho ds scale with the n umber of tags and the 13 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 SRAE JRAC 2D MUSIC 2D FFT IR First (a) Range est. error, 4-tag setup 0 10 20 30 40 50 60 0 0.2 0.4 0.6 0.8 1 SRAE JRAC 2D MUSIC 2D FFT (b) Angle est. error, 4-tag setup 0 2 4 6 8 10 Positioning Error (m) 0 0.2 0.4 0.6 0.8 1 CDF ML: SRAE IRLS: SRAE ML: JRAC IRLS: JRAC ML: 2D MUSIC IRLS: 2D MUSIC ML: 2D FFT IRLS: 2D FFT (c) Localization error, 4-tag setup 0 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 SRAE JRAC 2D MUSIC 2D FFT IR First (d) Range est. error, 100-tag setup 0 20 40 60 80 100 120 140 160 180 0 0.2 0.4 0.6 0.8 1 SRAE JRAC 2D MUSIC 2D FFT (e) Angle est. error, 100-tag setup 0 2 4 6 8 10 Positioning Error (m) 0 0.2 0.4 0.6 0.8 1 CDF ML: SRAE IRLS: SRAE ML: JRAC IRLS: JRAC ML: 2D MUSIC IRLS: 2D MUSIC ML: 2D FFT IRLS: 2D FFT (f) Localization error, 100-tag setup Fig. 6. Joint range and angle estimation, and lo calization error on exp erimen tal setup with (a-c) 4 tags, and (d-f ) 100 tags. (c,f ) A : B denotes the p ositioning metho d A is applied to the ranges and AoAs which are estimated by metho d B . ML metho d denotes ML-based Gradient Ascent with Line Search. n um ber of measurements p er tag, both of whic h are higher in the 100-tag dataset. In contrast, IRLS requires only the in v ersion of a tall matrix with xed rank, so its runtime gro ws m uch more slo wly . Ov erall, ML with gradient ascent demands signicantly more computational time: at least 10 × higher than IRLS in the 4-tag setup and 300 × higher in the 100-tag setup. These results highligh t that while ML and IRLS achiev e comparable p ositioning accuracy , IRLS oers a clear adv antage in scalability . Ov erall, even though the ranking of range-angle estima- tion algorithms diers b etw een range and angle accuracy , the ultimate metric of interest is positioning qualit y . Our results show that JRA C and SRAE, combined with IRLS, deliv er the b est trade-o b et w een accuracy and computa- tional eciency . This makes them strong candidates for practical m ulti-static lo calization systems. VI I. Conclusion In this w ork, w e introduced the JRAC and SRAE metho ds for real-time join t range–AoA estimation in m ul- tistatic BNs. These tw o metho ds achiev e range and angle measuremen t accuracy comparable to FFT- and subspace- based baselines, while signican tly reducing computa- tional complexity . Building on these measuremen ts, we in v estigated tag localization using b oth ML and IRLS algorithms. Our exp erimental results revealed that an unmo died IRLS approach, as commonly adopted in prior literature, can assign articially high weigh ts to outlier measurements, leading to degraded lo calization p erformance. T o address this limitation, we proposed a scaling sc heme that eectiv ely mitigates the impact of outliers. W e further demonstrated that the scaled IRLS algorithm ac hiev es localization accuracy comparable to brute-force searc h on the ML criterion, while signi- can tly reducing computational complexity . This reduction enables practical real-time tag lo calization, making the prop osed approach well suited for use in large-scale, time- critical BN deploymen ts. References [1] K. Zheng, R. Xu, J. Mei, H. Y ang, L. Lei, and X. W ang, “Ambien t IoT tow ard 6G: Standardization, potentials, and challenges,” IEEE Access, v ol. 12, pp. 146 668–146 677, 2024. [2] N. V an Huynh, D. T. Hoang, X. Lu, D. Niyato, P . W ang, and D. I. Kim, “Ambien t bac kscatter communications: A con tempo- rary surv. ” IEEE Commun. Surveys & T utorials, v ol. 20, no. 4, pp. 2889–2922, 2018. [3] M. M. Butt, N. R. Mangalvedhe, N. K. Pratas, J. Harreb ek, J. Kimionis, M. T ayy ab, O.-E. Barbu, R. Ratasuk, and B. V e- jlgaard, “Ambien t IoT: A missing link in 3GPP IoT devices landscape,” IEEE In ternet of Things Mag., vol. 7, no. 2, pp. 85–92, 2024. [4] Global, EPC, “EPC radio-frequency identit y proto cols class-1 generation-2 UHF RFID proto col for communications at 860 MHz–960 MHz,” GS1, 2008. [5] J. D. Grin and G. D. Durgin, “Complete Link Budgets for Backscatter-Radio and RFID Systems,” IEEE Antennas Propag. Mag., vol. 51, no. 2, pp. 11–25, 2009. [6] H. G. W ang, C. X. P ei, and C. H. Zhu, “A link analysis for passive UHF RFID system in LoS indoor environmen t,” Int. Conf. Wireless Commun., Net w. Mobile Comput., pp. 1–7, 2008. [7] J. Kimionis, A. Bletsas, and J. N. Sahalos, “Increased range bistatic scatter radio,” IEEE T rans. Commun., vol. 62, no. 3, pp. 1091–1104, 2014. [8] G. V ougioukas, S.-N. Daskalakis, and A. Bletsas, “Could battery-less scatter radio tags achiev e 270-meter range?” in IEEE wireless p ow er transfer Conf., 2016, pp. 1–3. [9] P . N. Alevizos, K. T ountas, and A. Bletsas, “Multistatic Scatter radio sensor netw orks for extended coverage,” IEEE T rans. Wireless Commun., v ol. 17, no. 7, pp. 4522–4535, 2018. [10] T. E. Abrudan, K. Patel, J. Kimionis, T. Esmaeilbeig, E. Kampianakis, S. D. Liyanaarac hc hi, and M. Eggleston, “Next-generation backscatter netw orks for integrated commu- nications and RF sensing,” arXiv preprint 2025. 14 [11] K. Patel, J. Zhang, J. Kimionis, L. Kampianakis, M. S. Eggleston, and J. Du, “Analyzing the scalability of bi-static backscatter netw orks for large scale applications,” IEEE J. of Radio F req. Identication, v ol. 9, pp. 6–16, 2025. [12] E. Kampianakis, J. Kimionis, K. T ountas, C. K onstan topoulos, E. Koutroulis, and A. Bletsas, “Wireless environmental sensor netw orking with analog scatter radio and timer principles,” IEEE Sensors J., vol. 14, no. 10, pp. 3365–3376, 2014. [13] G. V annucci, A. Bletsas, and D. Leigh, “A softw are-dened radio system for backscatter sensor netw orks,” IEEE T rans. Wireless Commun., vol. 7, no. 6, pp. 2170–2179, 2008. [14] D. Bharadia, K. R. Joshi, M. K otaru, and S. Katti, “Bac kFi: High throughput WiFi backscatter,” ACM SIGCOMM Comput. Commun. Rev., vol. 45, no. 4, pp. 283–296, 2015. [15] B. Kellogg, V. T alla, S. Gollakota, and J. R. Smith, “Passiv e Wi-Fi: Bringing low p o w er to Wi-Fi transmissions,” USENIX Symp. on Netw ork ed Syst. Design and Implementation, 2016. [16] P . Zhang, D. Bharadia, K. Joshi, and S. Katti, “HitchHike: Practical bac kscatter using commodity WiFi,” Pro c. A CM Conf. embedded netw ork sensor Syst. CD-ROM, pp. 259–271, 2016. [17] P . Zhang, C. Josephson, D. Bharadia, and S. Katti, “F reeRider: Backscatter communication using commo dity radios,” Proc. In t. Conf. Emerg. Netw. Experiments and T echnologies, pp. 389– 401, Nov. 2017. [18] B. Kellogg, A. Parks, S. Gollakota, J. R. Smith, and D. W ether- all, “Wi-Fi backscatter: Internet connectivity for RF-p ow ered devices,” A CM SIGCOMM Comput. Commun. Rev., vol. 44, no. 4, pp. 607–618, Aug. 2014. [19] V. T alla, M. Hessar, B. Kellogg, A. Na ja, J. R. Smith, and S. Gollakota, “LoRa backscatter: Enabling the vision of ubiq- uitous connectivity ,” Pro c. ACM Interactiv e, Mobile, W earable and Ubiquitous T echnologies, v ol. 1, no. 3, pp. 105:1–105:24, Sep. 2017. [20] S. W ang, Y. Y an, F. Han, Y. Tian, Y. Ding, P . Y ang, and X.-Y. Li, “MultiRider: Enabling multi-tag concurrent OFDM backscatter by taming in-band interference,” in Pro c. Annu. Int. Conf. Mobile Syst., Appl. and Services, 2024, pp. 292–303. [21] K. Skyv alakis, E. Giannelos, E. Andrianakis, and A. Bletsas, “Elliptical DoA estimation & lo calization,” IEEE J. of Radio F req. Iden tication, vol. 6, pp. 394–401, 2022. [22] M. V estakis, P . N. Alevizos, G. V ougioukas, and A. Bletsas, “Multistatic Narro wband Lo calization in Backscatter Sensor Netw orks,” in IEEE Int. W orkshop Signal Pro cess. Adv ances Wireless Commun., 2018, pp. 1–5. [23] J. Jiang, J. W ang, Y. Chen, S. T ong, P . Xie, Y. Liu, and Y. Liu, “Willow: Practical WiFi backscatter lo calization with parallel tags,” in Proc. Annu. In t. Conf. Mobile Syst., Appl. and Services, 2024, pp. 265–277. [24] J. Kim, J. Ch un, and S. Song, “Join t range and angle estimation for FMCW MIMO radar and its application,” arXiv preprint arXiv:1811.06715, 2018. [25] I. Jaafar, R. Amara, H. Boujemâa, and M. Siala, “Joint angle and delay estimation of p oint sources,” in Int. Conf. Electron. , Circuits and Syst. IEEE, 2005, pp. 1–5. [26] Y. Sc hröder and L. W olf, “InPhase: Phase-based ranging and localization,” ACM T rans. Sensor Net w., vol. 18, no. 2, 2022. [27] X. Xu, T. Huang, X. Kuai, and Y.-C. Liang, “Joint localization and signal detection for ambien t bac kscatter communication systems,” IEEE T rans. Wireless Commun., vol. 23, no. 10, pp. 14 437–14 451, 2024. [28] M. W u and C. Hao, “Super-Resolution T oA and A oA estimation for OFDM radar systems based on compressed sensing,” IEEE T rans. A erosp. and Electron. Syst., v ol. 58, no. 6, pp. 5730–5740, 2022. [29] M. Kotaru, K. Joshi, D. Bharadia, and S. Katti, “Sp otFi: Decimeter level lo calization using WiFi,” in Proc. ACM Conf. special interest group data Commun., 2015, pp. 269–282. [30] K. V. Mardia and P . E. Jupp, Directional Statistics. Wiley Series in Probability and Statistics, 1999. [31] T. E. Abrudan, Z. Xiao, A. Markham, and N. T rigoni, “Un- derground i ncremen tally deployed magneto-inductive 3-D posi- tioning net work,” IEEE T rans. Geosci. Remote Sens., vol. 54, no. 8, pp. 4376–4391, 2016. [32] S. W ang, B. R. Jackson, and R. Inkol, “Performance c har- acterization of AoA geolo cation systems using the von Mises distribution,” in IEEE V eh. T echnol. Conf., 2012, pp. 1–5. [33] H. Naseri and V. Koivunen, “A Bayesian algorithm for dis- tributed netw ork lo calization using distance and direction data,” IEEE T rans. Signal Inf. Pro cess. ov er Netw., vol. 5, no. 2, pp. 290–304, 2018. [34] K. P anw ar, P . Babu, and P . Stoica, “Maximum likelihoo d algorithm for time-delay based multistatic target localization,” IEEE Signal Pro cess. Lett., vol. 29, pp. 847–851, 2022. [35] F. Perez-Cruz, C.-K. Lin, and H. Huang, “BLADE: A universal, blind learning algorithm for T oA lo calization in NLoS channels,” in IEEE Glob ecom W orkshops, 2016, pp. 1–7. [36] C. Geng, T. E. Abrudan, V.-M. Kolmonen, and H. Huang, “Experimental study on probabilistic T oA and AoA joint lo- calization in real indoor environmen ts,” in IEEE Int. Conf. Commun., 2021, pp. 1–6. [37] Y. T. Chan and K. C. Ho, “A simple and ecient estimator for hyperb olic lo cation,” IEEE T rans. signal Pro cess., vol. 42, no. 8, pp. 1905–1915, 1994. [38] F. Ma, L. Y ang, M. Zhang, and F.-C. Guo, “TDoA source positioning in the presence of outliers,” IET Signal Pro cess., vol. 13, no. 7, pp. 679–688, 2019. [39] R. Amiri, F. Behnia, and H. Zamani, “Asymptotically ecien t target lo calization from bistatic range measuremen ts in dis- tributed MIMO radars,” IEEE Signal Pro cess. Lett., vol. 24, no. 3, pp. 299–303, 2017. [40] P . V arshney , P . Babu, and P . Stoica, “Outlier-robust m ultistatic target lo calization,” IEEE Signal Pro cess. Lett., 2025. [41] T. Abrudan, H. Claussen, and V.-M. Kolmonen, “Positioning system and metho d,” Nov. 19 2024, US Paten t 12,146,975. [42] O. Kanhere, S. Goyal, M. Beluri, and T. S. Rappap ort, “T arget localization using bistatic and multistatic radar with 5G NR wa veform,” IEEE V eh. T ec hnol. Conf., pp. 1–7, 2021. [43] N. Decarli and D . Dardari, “RFID and radar lo calization: A position error bound analysis,” in IEEE In t. Conf. Commun. W orkshops, 2013, pp. 37–41. [44] H. Marom, Y. Bar-Shalom, and B. Milgrom, “Bistatic Radar T racking With Signicantly Impro v ed Bistatic Range Accu- racy ,” IEEE T rans. A erosp. and Electron. Syst., vol. 59, no. 1, pp. 52–62, F eb. 2023. [45] T.-J. Shan, M. W ax, and T. Kailath, “On spatial smo othing for direction-of-arriv al estimation of coheren t signals,” IEEE T rans. Acoustics, Sp eech, and Signal Pro cess., vol. 33, no. 4, pp. 806– 811, 1985. [46] R. Schmidt, “Multiple emitter lo cation and signal parameter estimation,” IEEE T rans. An tennas and Propag., vol. 34, no. 3, pp. 276–280, 1986. [47] A. Bletsas, A. G. Dimitriou, and J. N. Sahalos, “Improving backscatter radio tag eciency ,” IEEE T rans. Microw. Theory and T echniques, vol. 58, no. 6, pp. 1502–1509, 2010. [48] J. Kimionis and M. M. T entzeris, “RF tag fron t-end design for uncompromised comm unication and harvesting,” in IEEE RFID T echnol. and Appl. Conf., 2014, pp. 109–114.