Markov Multi-Beamtracking on 60 GHz Mobile Channel Measurements

Beamforming training refers to the exhaustive scan over which the transmitter and receiver jointly steer their beams along a predeﬁned set of double-directional angles to determine the beam pairs that coincide with the dominant propagation paths of the channel, for spatial multiplexing at millimeter-wave. When mobile, training necessitates a high refresh rate to maintain connectivity and so, to reduce overhead, beamtracking algorithms exploit the spatial-temporal consistency of the channel to localize the scan around the beam pairs determined at a previous time. The algorithms’ true performance, however, is still unknown since results reported to date are based on oversimpliﬁed channel models. In this paper, we propose a novel beamtracking algorithm formulated as a ﬁrst-order Markov process that supports multiple beam pairs. The algorithm is evaluated through actual channel measurements – not a channel model – recorded with our high-precision 3D double-directional 60 GHz channel sounder. The measurement campaign, to our knowledge, is unprecedented: with 10, 895 large-scale measurements, spaced 8.8 cm apart on average to emulate continuous motion, over which the mobile receiver traversed a total of 900.2 m. We demonstrate that four beam pairs can be sustained always and that eight pairs can be sustained 57% of the time.


I. INTRODUCTION
The ever-rising demand for reliable and ubiquitous broadband access has prompted cellular providers to expand beyond the sub-6 GHz frequency bands -where 1G-4G networks have operated to date -to millimeter-wave (mmWave) bands for 5G -effectively 28-100 GHz -where 100x more bandwidth is available. The expansion, however, comes at the expense of greater free-space, penetration, and oxygenabsorption pathlosses. To compensate the link budget, highgain -and by corollary narrowbeam -antennas will be employed. Notably, the development of phased-array antennas (PAAs) at mmWave [1]- [3], which coherently phase the array antennas to synthesize pencilbeams -beams with array gain in excess of 18 dB and with beamwidth of just a few degrees -that are electronically steerable to provide omnidirectional coverage, has ushered in the 5G revolution.
The implication of shorter wavelengths is negligible diffraction due to the narrowing of the Fresnel zone [4], leaving the line-of-sight (LoS) path, and specular paths reflected from ambient scatterers as the dominant propagation paths in the mmWave channel. Accordingly, the transmitter (T) and receiver (R) beams -pairwise referred to as a doubledirectional (DD) beam -will be steered along the respective angle-of-departure (AoD) and angle-of-arrival (AoA) -pairwise referred to as a DD angle -of the dominant paths. To determine their DD angles, the channel must be estimated. For instance, the IEEE 802.11ay standard [5]- [7] for wireless local area network (WLAN) operating in the 60 GHz unlicensed band estimates the channel through a protocol known as beamforming training (BT): through dedicated pilot sequences, DD beams exhaustively scan a predefined set of DD angles, constituting a codebook, to identify the ones e.g. with the highest signal-to-noise ratio (SNR). The overhead associated with BT, however, is burdensome, diverting temporal and/or spectral radio resources that could otherwise be used for data transmission. While techniques based on compressed sensing can ease the overhead [8]- [12], they are applicable only in static or quasi-static channel conditions. What is more, in the presence of mobility, beams will quickly misalign with the dominant paths -a small angle misalignment with such narrow pencilbeams can inflict a huge drop in SNR -requiring a high refresh rate for BT, exacerbating the already burdensome overhead, especially in MIMO (Multiple-Input Multiple-Output) architectures, where multiple DD beams are supported through spatial multiplexing.
To deal with mobility, algorithms for beamtracking have emerged over the past five years [20]- [38]. The algorithms exploit the spatial-temporal consistency of the channel to reduce the extent of the BT scan. The underlying assumption is that the DD angle of the dominant paths will change incrementally with the incremental motion of the T, R, and/or ambient scatterers; so to realign the DD beams with the dominant paths while in motion, is it sufficient to scan angles within the codebook that are local to the DD beams determined at a previous time, effectively tracking the beams. Generally speaking, beamtracking algorithms exploit information about the channel from the past to localize the scan at present. In the sequel, we have categorized the main features of the algorithms cited above: r MIMO architecture: Most papers consider either analog beamforming [20]- [27] -where all antennas of a PAA are synthesized into a single DD beam to support a single digital chain in SISO (Single-Input Single-Output) -or hybrid beamforming -where subarrays of a PAA (or multiple PAAs) are synthesized into multiple beams to support multiple chains in SU (Single-User) or MU (Multiple-User)-MIMO [28]- [37]. Yet one paper [38] considers full digital beamforming, where each PAA antenna supports a single chain in MIMO, even though it bears a significantly higher implementation cost vs. hybrid beamforming. r number of hypotheses: Most papers commit to a single hypothesis per time [20]- [26], [28]- [32], [34], [35], meaning that only information about the DD beam that is estimated per chain is propagated over time, whereas some papers propagate information about multiple candidate DD beams per chain, either through particle filtering [27], [33], [36] or Bayesian inference [37], [38], and so are more robust to harsh, uncertain, or rapidly changing channel conditions. r tracking mechanism: Many papers assume an underlying motion model -typically linear -either directly for the DD angle of the dominant paths or vis-à-vis R's motion, and solve a least-squares fit between the angles predicted from the motion and the angles within the codebook, either through Kalman filtering [21], [23], [30], [34], gradient descent [24], [29], [31], or particle filtering [27], [33], [36]. Others, rather, localize the scan within the codebook by means of a probability distributionuniform [22], [28], Gaussian [24], [32], [37], or exponential [26], [38] -centered at a previously estimated angle; the codebook angles are then conditioned by the probability when generating a new estimate, e.g. through Bayesian inference [37], [38]. Yet other papers either simply widen/narrow the beamwidth [20], [25] or increase/decrease the BT refresh rate [35] based on R's speed.
r channel: Save one, all papers employ a channel model -not actual measurements. The most popular channel model is statistical, typically some variant of the spatial Saleh-Valenzuela model [21], [33] through which the channel is represented as a collection of paths clustered in the angle-delay domain per time; spatial consistency is imposed by correlating the channel over time through some motion model -linear, Brownian, etc. [23], [28], [29], [31], [35], [37], [38]. Other channel models are quasi-deterministic [27], [36], which combine a statistical model with raytracing such that spatial consistency is inherent, whereas other papers employ raytracing only [32]. Most papers assume a sparse channel, i.e. with just a couple of dominant paths [22], [30] or even just the LoS path [20], [24], [25], [26]. The single paper that employs actual measurements only considers motion over a meter or so and hence does not capture large-scale variation over which the number of scatterers and their properties change significantly. Indeed, the importance of accurate channel representation [39] cannot be underestimated for a meaningful evaluation of beamtracking algorithms. Although the LoS and specular paths account for most of the channel power, diffuse paths -which arise from rough surface scattering and so appear only when the wavelength is comparable to the dimension of the roughness, such as at mmWave -can account for up to 40% [40]- [44]. Diffuse paths tend to cluster around specular paths, forming dense pockets in angle (and delay). Consequently, when a beam is steered along a specular path, the surrounding diffuse paths temper the transition off the beam when in motion; at the same time, diffuse paths can act as a source of interference, not only within the beam, but also between different beams. Hence results generated from oversimplified models which assume that the channel can be accurately represented by just a few paths -when in reality, as we shall see, there can be hundreds -are unreliable. Moreover, all the cited papers assume a smooth trajectory of the specular paths as they traverse an object's surface, when in fact objects such as walls are composed from windows, doors, etc., which can have intricate multilayered profiles characterized by discontinuous surfaces made from composite materials (glass, metal, wood, etc.), giving rise to paths with non-specular properties that vary randomly in path gain, angle, and delay. Even the LoS path is not purely deterministic as it is subject to blockage. These realistic channel conditions can cause beams to easily lose track.
In this paper, we propose a beamtracking algorithm that combines the best features of the papers cited and evaluate its performance through actual measurements. Our three main contributions are as follows: 1) We recorded 10, 895 large-scale channel measurements with our high-precision 60 GHz switched-array channel sounder. The measurements were spaced 8.8 cm apart on average to emulate continuous motion, as our mobile receiver traversed a total distance of 900.2 m in a lecture room. Super-resolution techniques were then applied to extract distinct channel paths in path gain, delay, and 3D (in azimuth and elevation) DD angle from each measurement. 2) We formulate a beamtracking algorithm as a first-order Markov process that supports multiple beams through hybrid beamforming in SU-MIMO, that entertains multiple hypotheses, and that dynamically adjusts the scan locality within the codebook to the rate of motion. 3) We substantiate the robustness of the proposed algorithm by applying it to the channel measurements, demonstrating that up to four DD beams can be supported in the environment always, six beams 90% of the time, and eight beams 57% of the time. We also consider human presence. The remainder of the paper is composed as follows: in Section II, we describe our channel sounder and the extensive measurements campaign we conducted in the lecture room environment. Section III follows with our implementation of hybrid beamforming in SU-MIMO and then in Section IV with our proposed beamtracking algorithm. Section V reports the results of the algorithm applied to the channel measurements, and the last section is reserved for conclusions.

II. CHANNEL MEASUREMENTS
This section describes the channel measurements that were recorded for the purpose of evaluating our Markov multi-beamtracking algorithm. First, our channel sounder is presented, followed by details of the measurement campaign, then by the processing techniques implemented to extract channel paths and their properties from the measurements. Fig. 1(a) displays T and R of our 60 GHz 3D DD switchedarray channel sounder. R features a circular array of 16 horn antennas. The horns have a 3D Gaussian radiation pattern with 22.5°beamwidth and 18.1 dBi gain. To avoid "blind spots, " the angular spacing between adjacent horns on the array is matched to the beamwidth; specifically, the horns are offset by 22.5°in azimuth, and in elevation each other is offset by 22.5°; the resultant synthesized azimuth field-of-view (FoV) of the array is 360°and 45°in elevation. T is almost identical, except that it features a semicircular array of only 8 horns, limiting the azimuth FoV to 180°while maintaining the same 45°elevation FoV. Further details of the system are provided in [16].

A. CHANNEL SOUNDER
At T, an arbitrary waveform generator produces a repeating M-ary pseudorandom noise (PN) code, with 2047 chips and 2 GHz chip rate (0.5 ns delay resolution). The PN code is generated at IF, upconverted to an Radio Frequency (RF) center frequency of precisely 60.5 GHz, and then transmitted.
At R, the signal received is downconverted back to Intermediate Frequency (IF) and then sampled at 40 GHz (oversampled by 10x, both to reduce aliasing and to improve the SNR). Finally, the sampled signal is matched filtered with the PN code to generate a complex-valued channel impulse response (CIR) as a function of delay. For a T power of 100 mW, the maximum measurable path loss of the system is 162.2 dB when factoring in antenna gain, processing gain, system noise, and the other components of the link budget.
The PN code is electronically switched through each pair of T and R horns in sequence, resulting in 16 × 8 = 128 CIRs, which is referred to as an acquisition. Synchronous triggering between the T and R is provided through Rubidium clocks at each end, which are also used to discipline the local oscillators.

B. MEASUREMENT CAMPAIGN
The measurement campaign was conducted in a 19.3 m × 10 m lecture room, pictured in Fig. 1(a). T was fixed on a tripod at 2.5 m height and R was mounted on a mobile robot at 1.6 m height. The robot is equipped with an onboard computer that enables rapid and autonomous recording of channel acquisitions. It is also equipped with a laser-guided navigational system that reports its location and heading as follows: the robot first surveys the environment to generate a floor plan -pictured in Fig. 1(b) -through simulataneous localization and mapping (SLAM) [18] in a global coordinate system with centimeter accuracy. The location of the robot is reported in global coordinates while its heading is reported with respect to boresight of one of the 16 antennas on the R array and then transferred to global coordinates; that way the estimated DD angle of the paths is independent of the robot's heading. The top section of the room was occupied by two tables with surrounding chairs and the bottom section by two rows of chairs. The tables and chairs were wrapped with kraft paper to assist the robot with navigation.
During measurement, T was placed separately at the four corners of the room, marked T1-T4 in Fig. 1(b). Since its azimuth FoV is only 180°, T was oriented towards the center of the room, marked by an arrow. For each T location, R moved in a serpentine trajectory in each of the six areas, marked A-F, up to a pedestrian speed of 2 m/s while recording. The photograph in Fig. 1(a) was taken for Area F/T3. The average distance between consecutive recordings was 8.8 cm to emulate continuous motion. Between the four T locations and six R areas, a total of 10, 895 channel acquisitions was recorded. While recording, the channel was void of any motion (e.g. human) besides that of R, and LoS conditions were maintained throughout.

C. PATH EXTRACTION
The 128 CIRs per acquisition were coherently combined through the space-alternating generalized expectationmaximization (SAGE) super-resolution algorithm [17] to extract channel paths and their properties. Given the wide beamwidth of the horns (22.5°), the algorithm relies chiefly on the high delay resolution of our system (0.5 ns delay bin) to resolve different paths; notwithstanding, if two paths arrive in the same delay bin and within the same beamwidth, they can still be resolved to some extent if one of the paths is at least 6 dB stronger than the other (rule of thumb). Our SAGE algorithm is based on time difference of arrival (TDoA): If a path is detected in the same delay bin by a least three adjacent horns at T and by at least three adjacent horns at R, the relative difference between their delays against the relative displacement between their phase centers is used to estimate its DD angle. The SAGE algorithm is also based on power difference: the relative difference in the detected power of the path between the horns against the relative antenna gains along the path's DD angle is also incorporated to estimate its DD angle. Details of the algorithm are described in [19].
The output from SAGE for each acquisition was indexed against the recording time t, initialized to t = 0 s in each area. At each t, N (t ) channel paths indexed through n = 1 … N (t ) were extracted, complete with path properties in the six-dimensional domain: complex amplitude α n (t ), delay τ n (t ), and 3D DD angle denotes AoA, both in the global coordinate system of the environment map generated by the robot (see Fig. 1(b)). The measurement error of the channel sounder was computed against the ground-truth properties of the LoS path, whose AoD/AoA and delay were given from the geometry between T and R, and its path gain (|α n (t )| 2 ) from the delay through Friis transmission equation. The average measurement error over all areas and transmitters was reported as 1.1 dB in path gain, 0.54 ns in delay, and less than 3.7°in any angle dimension.
Any measurement taken with a channel sounder will contain not only the response of the channel, but also the response of the sounder itself, namely the directional patterns of the antennas and the responses of the T and R front ends. While the SAGE algorithm accounts for de-embedding the antenna patterns, the effects of the T and R front ends were removed through pre-distortion filters designed from a back-to-back calibration method [16]. Hence the path properties reflect the "pristine" response of the channel alone and not that of the measurement system. Fig. 2 displays the 349 paths measured at t = 0 s in Area D/T1, marked accordingly in Fig. 1(b). The paths are displayed in azimuth AoD vs. azimuth AoA vs. delay and are color-coded against path gain in the colorbar. As is evident from the plot, the paths cluster naturally in the angle-delay domain. The specular path -the strongest -of each cluster identified (through the algorithm in [41]) is marked by a square. The ambient scatterers that generated the specular paths were classified through inverse raytracing [42]; the scatterers classified are color-coded against the legend and are also marked in Fig. 1(a) and (b). Note that up to three reflections were classified (e.g. Left + Right + Left wall).
Although LoS conditions were maintained throughout while recording, the limited 180°azimuth FoV at T and the  Fig. 1(b). The paths are displayed in the azimuth AoD vs. azimuth AoA vs. delay and are color-coded against path gain in the color bar. It can be noted that the paths cluster naturally in the angle-delay domains. The specular path of each cluster identified is marked with a square that is color-coded in the legend against the scatterer that generated the cluster.
limited 45°elevation FoV at both T and R acted as to block the channel paths at times. Furthermore, the serpentine trajectory of the robot was chosen deliberately to generate a rapidly changing channel as the robot turned 180°in heading at the edges. When turning, the measurements error was up to twice as great than in the linear segments of the trajectory due to the robot's navigation system, which averages out heading errors over time, but as a result reports heading that lags the true heading -albeit by just a few degrees -yet significant enough to affect the precision of SAGE. Blockage, rapidly changing heading, and measurement error, together with an actual physical environment characterized by scatterers with intricate geometry and diverse materials, created channel conditions that rigorously tested the beamtracking algorithm.

III. HYBRID BEAMFORMING IN SU-MIMO
In full digital beamforming, each antenna supports a unique digital chain, meaning that the channel matrix is computed between all T antennas and all R antennas in the design of the precoder/combiner. While full digital beamforming can deliver optimal performance in terms of capacity, it is prohibitive when the number of antennas is massive, as is the case for mmWave PAAs, due to the messaging overhead to populate a massive channel matrix in addition to separate digitalto-analog converters/analog-to-digital converters (DAC/ADC) per chain converters per chain. To circumvent the problem, IEEE 802.11ay has provisioned for hybrid beamforminga combination of analog and digital beamforming -where each chain is associated with a subarray of PAA antennas (or a separate PAA). In analog beamforming, phase shifts are applied to the antennas to synthesize beams at T and R that are electronically steered along the DD angle of dominant channel paths -we henceforth refer to a T and R beam pair simply as a DD beam. Because the DD beams are narrow, the correlation between beams can be minimized -but not completely eliminated -through proper design. Digital beamforming then orthogonalizes any residual correlation by precoding the beams at T and combining the beams at R.
Our Markov multi-beamtracking algorithm considers SU-MIMO to support M DD beams. Given the reduced channelstate information compared to digital beamforming, hybrid beamforming will in general be suboptimal. Many hybridbeamforming algorithms have emerged in recent years - [15] provides a comprehensive survey -and performance will greatly depend on the channel. In this section, we propose a hybrid beamforming scheme that integrates our channel measurements directly, with no need for a channel model.

A. ANALOG BEAMFORMING
BT is realized through analog beamforming. Specifically, DD beams synthesized at T and R scan a set of predefined DD angles [θ T i θ R j ] composed from the Cartesian product of L x L single-directional angles, θ T i , i = 1 . . . L and θ R j , j = 1 . . . L respectively, with uniform spacing θ in between at each end, composing a codebook of size L 2 . Since the horns we used for measurement are not electronically steerable, in this paper we employ a planar PAA model instead, represented by the P-length steering vector s(θ) = 1 √ P e j 2π λ (cos θ A cos θ E ,sin θ A cos θ E ,sin θ E )·X , where X is 3 × P matrix of the 3D locations of the P antennas [45]. The beamformed CIR per codebook entry i j is then given from the measured at time t as [5]: where † denotes the Hermitian.

B. DIGITAL BEAMFORMING
Once the L 2 DD beams within the codebook have been scanned, the next step is to find the combination of M DD beams that optimizes SU-MIMO performance. To that end, we consider the K = L 2M possible permutations of the DD beams, permuted into M-dimensional sets indexed through Each set k has a digital channel matrix composed from the beamformed CIRs between the M x M single-directional beams in the set: (2) Hence the channel matrix also accounts for the interference between the M DD beams in the set. In fact, the metric we use to quantify the optimality of a set is through its ergodic channel capacity where ρ/M is the total T power set to 1 mW and σ 2 is the noise power at R, equivalent to 8 × 10 −9 mW for the 2 GHz bandwidth considered. Essentially, the optimal set is the one with the highest signal-to-interference-plus-noise ratio (SINR), i.e. the one that minimizes interference between the M DD beams. As such, capacity is a metric that considers the beams dependently.

IV. MARKOV MULTI-BEAMTRACKING
Our Markov multi-beamtracking algorithm is initialized through BT, through which all L 2 DD beams within the codebook are scanned exhaustively and the set of M permuted DD beams that yields the highest capacity is selected for hybrid beamforming. In principle, BT could be performed at each time, but, aside from the cumbersome overhead, the scanned beams are subject to channel uncertainty -rapid channel variation due to small-scale fading or short-time shadowing in cluttered environments, channel estimation error, etc. -compromising the stability of the beams over time. Unstable beams in turn trigger more frequent BT, forming a vicious circle. Our objective, rather, is to minimize the amount of channel scanning to reduce overhead while fusing the scans over time to foster beam stability. How we achieve this objective -inspired by an algorithm from robot localization and tracking [13] -is described in this section.

A. FIRST-ORDER MARKOV PROCESS
We pose our multi-beamtracking algorithm as a discrete firstorder Markov process. In general, a Markov process is defined by a discrete number of states, state observations over time, and fixed state transition probabilities between two consecutive times. In our algorithm, the states are the sets k , k = 1 . . . K and the observations are the capacity of the sets . The Markov process is then governed through the following equation [13]: The prior probability p(C(t )| k ) represents the likelihood of set k given the observation at a single time alone (at t); the set is represented as a probability to reflect the channel uncertainty in the observation. To favor sets with greater capacity, they are rewarded higher probabilities. Specifically, the prior probability of k is scaled according to its capacity C k (t ) with respect to the total capacity over all K sets as: Naturally, the prior probabilities over all K sets sum to 1. The posterior probability p( k |C(0) . . . C(t − 1)) represents the likelihood of set k given the observations over all times, from t = 0 s, when BT occurs, up to the previous observation at time t − 1. When a new capacity C(t ) is available at t, the posterior probabilities of all K sets are propagated to from t − 1 to t through (4). The term normalizes the posterior probabilities after propagation, imposing the law of total probability such that the sum over all sets is 1. The prior and posterior probabilities work hand-in-hand: the former admits changing channel conditions while the latter resists change, fostering beam stability over time. And as opposed to most other greedy algorithms [20]- [26], [28]- [32], [34], [35], the Markov process supports multiple hypotheses: it propagates all K states over time instead of committing to only the most likely. This allows hypotheses that persist over time to emerge, in which the prior reinforces the posterior, rendering the algorithm robust. Finally, because the capacity metric accounts for the interference between the M DD beams, the beams are tracked dependently, as opposed to the greedy alternative in which the beams are tracked independently of each other [30], [31]. Fig. 3 displays the measured paths (in azimuth AoD vs. azimuth AoA only) in Area D/T1 at seven times, marked t = 0:25:150 s in Fig. 1(b), over which the robot moves down the room, turns, and then moves back up, forming a "U" trajectory. The side-by-side subfigures t = (0, 150) s, t = (25, 125) s, t = (50, 100) s correspond to the side-by-side times on the "U", where the robot is roughly at the same location; in fact, note that the paths at the side-by-side times look very similar. In particular, note that the paths change DD angle incrementally over time and, more specifically, that all clustered paths move in unison. Because the channel is highly correlated in angle over time, it is expected that any DD beam will change DD angle incrementally between two consecutive times.
This spatial-temporal correlation of the channel is incorporated in the Markov process through the transition probabilities p( k | k ), to foster beam stability against channel uncertainty. Namely, set k can transition to k only if k is local to k as defined by: where r is an integer radius parameter demarcating the scan locality as multiples of the codebook spacing θ. Essentially, any set k can transition to set k only if k lies within the 3D DD (four-dimensional) cuboid centered at [θ T km θ R km ] = [θ T,A km θ T,E km θ R,A km θ R,E km ] per beam m. The number of beams in the cuboid is (3r) 4 . Fig. 4 illustrates the transition probabilities. As we shall see later, the radius r dynamically adjusts to how rapidly the channel changes over time, to reduce overhead.
The set MLE with the highest posterior probability is the Maximum Likelihood Estimation (MLE) used for hybrid beamforming at t: The MLE set is found through an exhaustive search. How to reduce the search space for practical implementation is discussed next.

B. IMPLEMENTATION
The theoretical aspects of our tracking algorithm have been described thus far; the practical aspects of the implementation are described in the sequel. In this paper, each beam is synthesized by an 8 x8 PAA subarray at both T and R with half-wavelength spacing between antennas (λ = 4.96 mm at 60.5 GHz), and equivalent array gain of 18 dB. We also assume θ = (12°, 12°), which is about the half-power beamwidth of the subarray, requiring L = 360 • 12 • · 360 • 12 • = 900 single-directional scans, or L 2 = 81, 000 DD scans to provide 3D omnidirectional coverage.

1) BEAMFORMING TRAINING (BT)
In theory, the L 2 DD beams scanned are permuted into K = L 2M sets of beams; in practice, K be a very large number. To reduce the number of sets permuted, we apply an interference constraint: no two beams in a set can have the same AoD nor AoA; their beams would otherwise interfere with each other, resulting in a set with low SINR (and in turn low capacity) anyway. While the interference constraint reduces the number of permutations down to K = ( L! (L−M )! ) 2 , this number in general will still be too large to permute. To reduce the number further, the L 2 beams are first filtered by SNR. The filtering exploits the sparsity of the mmWave channel. Namely, as illustrated in Fig. 3, the measured paths form narrow, densely packed clusters in the DD angle domain, therefore relatively few beams will actually encounter paths during BT, meaning that most beams will have negligible SNR. As an example,

2) BEAMTRACKING
After BT, the posterior probabilities are pinned to the prior probabilitiesi.e. p( k |C(0)) ≡ p(C(0)| k ) -given that only observations at a single time (t = 0 s) are available. Accordingly, the initial posterior probabilities are computed through (5) and propagated forward through (4). While the SNR filter and the interference constraint limit the number of initialized sets, that number is multiplied during propagation via the sum in (4), meaning that setk at t − 1 can transition to any other set k at t so long as the transition probability between them, p( k | k ), is non-zero. Thanks to our binary transition probabilities, the propagation is limited to the sets k that lie within the M cuboids of setk. Nevertheless, the number of sets can grow exponentially through propagation, by up to M · (3r) 4 per time. In order to curb growth, the number of sets propagated forward is limited to K MAX based on the CONDENSATION technique in [14]. In the technique, the K MAX sets with the highest posterior probabilities (out of the maximum K MAX · M · (3r) 4 propagated) are kept while the others are discarded; accordingly, the normalization term η(t ) in (6) is based on K MAX . Fig. 6 shows the 20 sets with the highest posterior probabilities over time in Area D/T1 for M MAX = 1; the value  Fig. 1(b), together forming a "U" trajectory by the robot motion. The paths are displayed in azimuth AoD vs. azimuth AoA and are color-coded against path gain in the colorbar. The four DD beams determined by the beamtracking algorithm are displayed as diamonds, color-coded against the scatterers classified in the legend, against which the beams are steered. As the robot moves down the "U", the beams move in lockstep with the DD angle of the scatterers. At t = 75 s, when the robot turns, the algorithm loses track of the Glass door scatterer and subsequently of the Right wall at 100 s. On the way back up the "U", the algorithm converges to a different set of scatterers, without triggering BT thanks to the multiple hypotheses.  M MAX = 1 was chosen deliberately to illustrate the dynamics of the single beam. At t = 0 s, the robot is farthest from T1 and although the LoS path is the strongest, the single beam latches on to the Glass Door scatterer instead, which generates rich diffuse paths that, when combined with the specular path, give rise to sets (DD beams) with the greatest capacities (SNR) by far; in fact, the posterior probabilities are concentrated in just those few sets. As the robot moves down the "U" and, in particular at t = 75 s, the reflection from the Glass door goes undetected -it is no longer "visible" due to lack of incidence at those T-R locations -so the remaining scatterers are much more comparable in strength; in turn, the posterior probabilities are spread out. The cyclical pattern in the posterior probabilities corresponds to the three consecutive "U"s the robot traverses in the area. In this paper, K MAX was determined as the number of sets with posterior probability within 99% of the highest.

3) DYNAMIC SCAN LOCALITY
As explained earlier, the binary transition probabilities limit the number of sets propagated over time. This, in addition to K MAX , eases the computational burden of maintaining potentially millions of sets. But where the limitation is most critical is in reducing overhead, by localizing the scan within the codebook from L 2 DD beams to just K MAX · M · (3r) 4 DD beams, freeing up radio resources for data transmission. The scan locality is adjusted dynamically through r based on how quickly the channel is changing: After BT, the cuboid radius is initialized to r = 1. A sharp drop in the MLE capacity indicates that at least one of the DD beams has fallen out of its own cuboid. Accordingly, if the MLE capacity at time t drops below a certain percentage with respect to t − 1 -in this paper we use 80%r is incremented by 1 per time until all the beams lie within the expanded cuboids, indicated by when the capacity rises back above the 80% baseline. The increment is limited to r ≤ 3, beyond which the multi-beamtracking is most likely irrecoverably lost; in that case, BT is triggered. If, however, the baseline capacity can be restored, r is decremented by 1 per time so long as the baseline capacity is maintained, until the equilibrium radius r = 1 is reached. Fig. 7 illustrates the dynamic radius over time in Area F/T4 for M MAX = 4: any steep drop in capacity triggers an expansion from r = 1 -at times up to r = 3 -and an eventual contraction back to equilibrium. The mechanism was successful in tracking all four beams throughout the whole area, without triggering BT.

V. RESULTS
In this section, we present the results of applying our Markov multi-beamtracking algorithm to the measurements recorded. The algorithm is run separately on each area/transmitter since the measurements were recorded disjointly. The first three subsections analyze performance at select times in an area, over all times in the area, and over all areas/transmitters, respectively. The final subsection, rather, analyzes performance in the presence of humans over all areas/transmitters, by applying a blockage model to the measurements.

A. SELECT TIMES IN AREA D/T1
Consider once more the measured paths in Fig. 3 for times t = 0:25:150 s in Area D/T1, over which the robot moves in the "U" trajectory shown in Fig. 1(b). Superimposed on the paths are the four DD beams determined by the tracking algorithm for M MAX = 4, displayed as diamonds color-coded against the scatterers classified in the legend, along which the beams are steered. Although, as noted before, the paths between the side-by-side times t = (0, 150) s, t = (25,125) s, and t = (50,100) s -at which the robot is roughly at the same location -are very similar, the DD beams tracked at the side-by-side times are not the same: From t = 0-50 s, the four beams reliably track the LoS path, Ceiling, Glass door, and Right wall scatterers. At t = 75 s, however, the Glass door is no longer visible, as explained earlier, and so its beam loses track. Yet, by virtue of multiple hypotheses, the algorithm latches on to the Right wall + Glass door instead. As the robot moves back up the "U" from t = 100-150 s, the Glass door is reabsorbed, but the Right wall is eventually replaced by the Left wall; in fact, by t = 150 s the purple diamond is perfectly aligned with the Left wall cluster. All the while, four beams are tracked over the whole "U" -even around the turn at t = 75 s when the channel is changing rapidly and the angle measurement error is twice as great as the rest -thanks to multiple hypotheses and dynamic scan locality.

B. ALL TIMES IN AREA D/T1
We now consider the robot motion over the whole trajectory in Area D/T1. Of course, as displayed in Fig. 8(c), the MLE capacity is also progressively greater from M MAX = 4 to 8 thanks to more and more beams. It is interesting to observe what happens when the number of sustainable beams across M MAX = 4, 6, 8 is equal; this happens only when M = 4, when the robot turns on the "U". There the algorithm finds it difficult to sustain more than four beams due to the rapidly changing channel and larger measurement error. In fact, in the attempt to support more beams (for M MAX = 6, 8), the algorithm triggers BT repeatedly but actually obtains worse capacity than when the algorithm is tracking (for M MAX = 4). This means that the suboptimal solution provided through BT with SNR filtering (for M MAX = 6, 8) -optimality can only be guaranteed without SNR filtering, but then the computational expense is prohibitive -is actual worse than the suboptimal solution provided through beamtracking with local scanning (for M MAX = 4).

C. ALL AREAS/TRANSMITTERS
The dynamics observed in Area D/T1 in fact were found to apply everywhere, as is apparent from the histograms in Fig. 9, computed over all areas and transmitters. Fig. 9(a) displays the distribution in the number of sustainable beams: for M MAX = 4 four beams are always sustainable, while for M MAX = 6 six beams are sustainable 90% of the time, and for M MAX = 8 eight beams are sustainable only 57% of the time. Fig. 9(b) and (c) respectively display the average overhead and average MLE capacity per M: for M MAX = 6, 8, the overhead increases with M due to the broader scan locality inherent to more beams as well as due to the fact that the additional beams are weaker and so lose track more easily, triggering BT more often; however, the capacity does increase with more beams, albeit at diminishing returns since the additional beams are weaker. To compare across M MAX = 4, 6, 8, the mean of the average overhead and average capacity over all M (weighted by the occurrence of each) are displayed in Fig. 9(b) and (c) as dashed lines. Again, the capacity increases with M MAX , but at the price of increasing overhead.

D. HUMAN PRESENCE
To analyze performance in the presence of humans, we applied the Modified Double Knife Edge Diffraction (MDKED) model described in [46], which is a blockage model reduced from measurements we conducted on human subjects using the same 60 GHz channel sounder. In the model, the human is represented as an infinitesimally high cylinder with 40 cm diameter; in our application, the cylinder is placed randomly in the lecture room, and based on the placement will block some of the measured paths. If a measured path is blocked, the path is diffracted around the cylinder, which is equivalent to adding an attenuation predicted by the blockage model to the measured path gain. As an example, Fig. 10(a) shows the measured path gain of the LoS path as the robot traverses Area D/T1, with and without human blockage. In the latter case, the measured path gain follows the periodic trajectory of the robot, while in the former the measured path gain is attenuated as the LoS path is gradually blocked by the cylinder -reaching a trough when the path intersects its center (see inset) -and eventually retreats back to equilibrium.
We consider three cases in our analysis, for H = 5, 10, 15 humans placed randomly and for M MAX = 4 fixed. The overhead and capacity averaged over all the areas/transmitters are shown in Fig. 10(b) and (c), respectively. Note that M = 4 beams can still be sustained always -even for H = 15and that when compared to the no human case in Fig 9(b) and (c), overhead increases while capacity decreases, as expected. Unexpected is that as H increases, the overhead actually decreases. This is because, when the room is densely packed (H = 15), many of the beams are blocked simultaneously and so propagation by way of numerous diffracted paths is prevalent, resulting in a posterior probability that is spread over multiple sets; the benefit to overhead is that the multiple sets propagate readily over time, seldomly triggering expansion of the cuboid radius or BT. Rather, when the human presence is light (H = 5), few beams are blocked simultaneously and so the posterior probability converges to just a few sets; when blockage does occur, the radius must expand significantly to maintain track at the onset of severe blockage, or BT is triggered. Of course, propagation via diffraction is weaker, in turn decreasing the SINR and overall capacity for H = 15 vs. H = 5.

VI. CONCLUSION
Beamtracking algorithms have emerged over the last five years to reduce the overhead inherent to the radio-resourceintensive beamscanning protocol that is employed for estimating the highly directional mmWave channel. Of the 19 papers on the topic that we found in literature, all but one employs an oversimplified channel model to gauge the performance of the algorithms, and the one paper that employs actual measurements considers small-scale motion only, over about a meter. The Markov multi-beamtracking algorithm that we propose in this paper, instead, was evaluated on a total of 10, 895 largescale channel measurements recorded with our high-precision 3D double-directional 60 GHz channel sounder, over which the receiver traversed a total of 900.2 m. The measurements were subject to a blockage, a rapidly changing channel, and measurement error together with an actual physical environment characterized by scatterers with intricate geometry and diverse materials, creating channel conditions that rigorously tested the proposed algorithm. Notwithstanding, the algorithm could sustain four beams always at an average capacity of 10.7 bit/s/Hz with an average overhead of 613 scans, six beams 90% of the time at 13.38 bit/s/Hz with 1317 scans on average, and eight beams 57% of the time at 20.2 bit/s/Hz with 2872 scans on average. In the presence of up to 15 humans in the environment, four beams could still be maintained, but capacity dropped 8.6 bit/s/Hz while overhead rose to 777 scans on average.