Outdoor Navigation Using Bluetooth Angle-of-Arrival Measurements

The introduction of direction finding in the Bluetooth standard enabled the use of antenna arrays for locating Bluetooth devices, using carrier phase measurements to estimate the direction from the array to a moving device. In this work, this feature is utilized for outdoor localization. We show how using repeated measurements from all array elements, instead of only the initial single-element reference samples as often suggested, can contribute to an improved estimate of the signal’s unknown carrier frequency offset, thereby improving the direction estimation performance. To run the direction-of-arrival estimation in real-time with high angular resolution on an embedded computer we propose a pseudo-spectrum peak search strategy that combines a coarse search, where the resolution is decided based on the array’s main lobe width, with a local nonlinear optimization for estimate refinement. We consider practical aspects relating to the phase sampling configuration and demonstrate direction estimation at up to 700m range with insignificant packet loss within 500m, and without significant loss of angular precision even when the received signal is near the receiver’s signal strength sensitivity threshold. In an open outdoor environment, using a square antenna array with 12 elements, the azimuthal performance is found to be very consistent with range, with noise standard deviation typically around 1°. While the elevation is significantly affected by multipath at lower elevation angles, with visible disagreement between frequency channels, it is shown to be consistent with simulations of ground reflection multipath.


I. INTRODUCTION
In 2019, the Bluetooth Special Interest Group (SIG) presented the Bluetooth 5.1 specification, which introduced support for direction finding using antenna arrays. This can provide estimates of the direction from the array to a moving tag in the antenna coordinate frame, but does not at the present time provide ranging. While the suggested applications typically focus on indoor tracking of devices, it can also be used outdoors where signal conditions are typically better, with fewer reflections off walls and other objects. Global Navigation Satellite System (GNSS) receivers are the most The associate editor coordinating the review of this manuscript and approving it for publication was Mohammad Zia Ur Rahman . common position source used for outdoor navigation, but the very weak signals are susceptible to interference such as jamming and ionospheric scintillation. As a result, having access to supplements or alternatives to GNSS is beneficial to ensure navigation robustness.
In 2019, [1] implemented the Bluetooth 5.1 direction finding using Software Defined Radios (SDRs) with a twoantenna array, since commercial off-the-shelf (COTS) equipment was not yet available, testing both outdoors and indoors. It was found that if the signal propagation direction is close to the linear array direction, phase delays appeared almost random. Different frequency channels performed differently, and averaging between them was found beneficial. We note that each packet is sent on a single channel, and that to average over several channels many packets must be used, which is a disadvantage since measurements on different channels will not be simultaneous.
In [2], COTS direction finding equipment from Texas Instruments was tested indoors and outdoors at short range, finding good performance in general, but significantly better outdoors. Two separate receiver arrays were used for 2D position estimation. Indoor experiments have also been conducted by others [3], [4], [5].
It should be noted that all previous work using Bluetooth 5.1 direction finding we have found focused on short-range indoor positioning applications. Outdoor testing was only performed for comparison with an environment offering cleaner signal propagation conditions, mainly with less multipath impact. This is very different from the outdoor applications we consider, demonstrating new capabilities for this low cost-equipment. Most previous experimental work has also focused on the estimation of a single azimuth angle for positioning in the horizontal plane using linear arrays or pairs of linear arrays such as the Texas Instruments BOOSTXL-AOA [2], [4], [6], [7], [8]. An exception is [3] where both a uniform rectangular array and an array of the same type as used in this paper were tested.
Outdoor navigation using phased array radio has been demonstrated in [9] using a system from Radionor Communications. This is a more expensive system primarily designed for long-range communication, with a range of tens of kilometers or more. The Radionor system estimates direction in order to perform active beamforming between the units, while Bluetooth direction finding only uses a single transceiver and hardware radio frequency (RF) signal switches on the array. Bluetooth is aimed at mass-market use, while the Radionor system uses licensed frequency bands. Similar to the Radionor system, Bluetooth can be used as a combined communication and navigation system. The Radionor system supports ranging, while this is not yet a part of the Bluetooth specification (see [10] for information about a planned ranging feature), although custom ranging solutions compatible with the standard exist, see e.g. [11].
Indoor short-range comparison of Ultra-wideband (UWB) and Bluetooth arrays for direction finding was performed in [7], finding that UWB in general had better performance, especially for complex environments with multipath and signal obstruction. The test did however use very small arrays with 2 or 3 elements, finding that 3 elements performed significantly better than 2 elements, suggesting that larger arrays could improve performance further. UWB uses a high sampling rate and one receiver per antenna, all with a common clock, and is therefore more expensive. The cost of additional array elements with Bluetooth is low, but as the number of elements increases the maximum number of measurements from each element is reduced due to the maximum length of the received signal. For UWB it is also common to position using trilateration with range measurements to multiple separate fixed anchors [12], [13], but this has the disadvantage of a requiring a more elaborate equipment setup than one or two arrays.
Bluetooth is typically considered a short-range system for use under 100 m distance, but a range of hundreds of meters is possible outdoors, and even kilometer level within the Bluetooth specification limits of transmit power, depending on the directionality of the antennas used. Bluetooth supports a mode offering increased range at the expense of reduced data rate, known as ''LE Coded PHY'', but this is not compatible with direction finding [14, p. 292]. The Bluetooth standard does not specify the direction estimation method, and many algorithms can be used for measurement processing. Examples of outdoor applications where Bluetooth direction finding could be used, if range and direction estimation performance is sufficient, are automatic docking of boats and ferries, robotic lawn mowers, and precision landing of unmanned aerial vehicles (UAVs).
The main contribution of this paper is the demonstration of long range use of Bluetooth direction finding outdoors in field experiments. Methods for • high precision carrier frequency offset estimation utilizing all available measurements • efficient high resolution spatial pseudo-spectrum peak search using optimization are proposed, resulting in improved estimation accuracy and precision of the direction of the received signal, and reduced computational load and latency. While the optimizationbased peak-search has not been used in the field experiments, it has been used in offline estimation for simulated multipath measurements, where the direction has been estimated for a large number of measurement sets with very high resolution, showing its usefulness.
The paper proceeds as follows: Section II explains the mathematical notation used throughout the paper. Section III introduces the Bluetooth direction finding feature, explains the measurements sampled by the receiver and how these can be used for Angle-of-Arrival estimation. Methods for high precision frequency deviation estimation and efficient high resolution pseudo-spectrum peak search are proposed. In Section IV, field experiments using a multirotor UAV and a fixed-wing UAV are conducted, demonstrating direction estimation out to 700 m range. Section V investigates the effect of ground reflection multipath interference on the direction estimate for vehicle navigation, and how different parameters affect the estimation error. Section VI concludes the paper with suggestions for future work.

II. PRELIMINARIES
The Euclidean norm is denoted · 2 with · being the variable placeholder.  The magnitude of a complex number is denoted | · |, the matrix transpose (·) , and the complex conjugate transpose (·) H . The imaginary unit is denoted j and the set of real and complex numbers R and C, respectively. v = v 1 v 2 is the element-wise product of the vectors v 1 and v 2 .

B. COORDINATES
The antenna array coordinate frame {a} is illustrated in Fig. 1. The navigation coordinate frame {n} has its origin coincident with {a}, but with axes pointing towards North-East-Down (NED). If the array is placed flat on the ground with x a towards North, the frames are identical. Directions in the antenna frame are parameterized using the polar angle α and the azimuthal angle . α is the angle of incidence, which is 0 in the boresight direction and (π/2) for a direction in the array xy-plane. is measured in the antenna xy-plane about z a using the right-hand rule, with = 0 for the direction x a . In {n} we have the azimuth angle n measured relative to North and elevation angle α n measured from the horizontal tangent plane.

III. BLUETOOTH DIRECTION FINDING
Bluetooth is a technical standard for short-range wireless communication using the 2.4-2.5 GHz industrial, scientific and medical (ISM) band. Bluetooth Low Energy (BLE), which is a subset of the Bluetooth standard intended for low power consumption, uses a total of 40 channels, from 2402 MHz to 2480 MHz, with 2 MHz channel spacing. The three channels at 2402 MHz, 2426 MHz, and 2480 MHz are used for advertising, while the rest are used for data transmission. Gaussian frequency-shift keying (GFSK) is used for modulation [10], which means that the frequency deviates with a positive offset for a binary 1, and a negative offset for a binary 0, but that a Gaussian filter is used on the baseband data before modulation to smooth the transitions and reduce sideband levels. For BLE using 1 Mb/s the nominal deviation is 250 kHz, and the average deviation should be between 225kHz and 275kHz, with 99.9% being above 185kHz [15].
The Bluetooth 5.1 Core specification introduced direction finding capability using antenna arrays. This is done by appending a Constant Tone Extension (CTE) to the Bluetooth packet transmitted, which is essentially a stream of binary ones, resulting in a sine wave at a fixed frequency at the end of the message. Receivers can measure the phase of this signal and perform carrier phase interferometry calculations. While antenna arrays are used either on the receive or transmit side, active beamforming is not used, and only a single antenna of the array is used at once. This allows low-cost hardware, using single-channel transceivers combined with electronically controlled RF switches. The sequential use of each antenna in the array does however mean that processing software is more complex compared to systems where all measurements are simultaneous.
Direction finding can use either Angle-of-Arrival, AoA, where the moving tag transmits the CTE and the array receiver measures phase, or Angle-of-Departure, AoD, where CTE transmission is switched between all array antennas, and the moving tag measures phase [16]. These two methods each have benefits and drawbacks. A disadvantage of AoD is that antenna switching occurs at the transmitting side, where significant transmit power must be switched. If multiple moving tags should be located, AoA requires all tracked devices to transmit data to the receiving array, and they should not transmit at the same time. Randomization of the transmit time can make the probability of traffic collision less likely, but with more devices, it can still occur. AoD only requires that the array transmits and all devices receive, allowing an unlimited amount of devices to estimate their direction from the array. If multiple arrays are to be used to track a single moving tag, however, AoA only requires the tag to transmit, while all arrays sample simultaneously. For AoD all arrays would need to transmit separately, meaning that direction estimates for each array would not be simultaneous. In this paper, only AoA is used, as it is the only feature supported by the equipment setup we are using for experiments. In the following, the terms moving tag and transmitter are therefore interchangeable.

A. CTE TRANSMISSION AND MEASUREMENT SAMPLING
The Bluetooth specification sets the maximum length of the CTE at 160 µs. The first 4 µs is the guard period, where no sampling is performed. This is followed by a 8 µs reference period, where one of the array elements is sampled repeatedly 8 times with 1 µs spacing. After the reference period, the remainder of the CTE is used to sample the array elements by using the signal switches to sequentially connect each element to the receiver. This period is divided into alternating switching slots and sampling slots of 1 µs or 2 µs duration [17]. Since a transient occurs when the receiver switches between antennas, the switching slots are intended to allow the signal to settle before sampling begins.
The transmitter broadcasts the CTE with a frequency f t = f c + f m + f t , where f c is the channel center frequency. The frequency deviation f m used for data modulation, nominally 250 kHz, is added since the CTE consists solely of digital ones. f t is the error in the channel center frequency due to the inaccuracy of the transmitter oscillator. The Bluetooth modulation frequency is not required to be very accurate, or constant over time, but can reasonably be assumed stable for the short duration of a single CTE. The receiver demodulates the signal using a local reference frequency set to the channel center frequency according to its own clock, f r = f c + f r , where f r is the receiver frequency error. For advertising channel 39 with a nominal frequency of 2480 MHz, where the CTE will have an intended frequency of 2480.25 MHz, the in-phase and quadrature (IQ) samples are generated using a local reference with a nominal frequency of 2480 MHz. This means that the IQ samples are not at baseband, but at an intermediate frequency offset from the baseband by approximately 250 kHz. This is known as the carrier frequency offset (CFO). The frequency deviation leads to a change in the phase relative to the receiver reference of approximately 90 degrees for each microsecond between samples. This is visible for the initial reference samples in Fig. 2.
Unlike receivers with a separate channel for each antenna in the array, allowing simultaneous sampling of all elements, BLE direction finding uses a single receiver connected to the array antennas using electronically controlled RF switches. The receiver samples all antennas sequentially by switching between them with known order and timing. Since measurements are not simultaneous and not at baseband, the frequency deviation must be estimated and taken into account in the direction estimation. Fig. 3 shows an example of samples from a CTE where the nominal 250 kHz has been used for compensation. There is still significant phase drift remaining, indicating that the actual deviation frequency is not 250 kHz. Compensating with an estimated 257 kHz gives the result in Fig. 4, where phases are close to constant. The carrier wavelength used can also be calculated from the estimated signal wavelength, but the error in this caused by a few tens of kilohertz error is very small. Because of the measurement noise present in the IQ samples, increasing the number of samples from a CTE used for frequency estimation can reduce the noise in the resulting frequency estimates, as will be explained in Section III-B3.

B. DIRECTION ESTIMATION FROM CTE IQ SAMPLES
Using the Bluetooth IQ measurements, the direction to the transmitter is estimated using the following steps: 1) Estimate the CFO of the signal using the classical periodogram [18] generalized to non-uniform sampling,  and in our case, multiple antennas. This creates a one-dimensional frequency spectrum. A frequency spectrum peak search is then used to determine the signal frequency. 2) Use the CFO estimate to correct the IQ measurements such that they can be treated as being simultaneous. 3) Create a spatial pseudo-spectrum using a conventional beamformer method [19], [20] with a steering vector based on the far-field assumption. 4) Find the direction estimate using a peak search of the spatial pseudo-spectrum.

1) STEERING VECTOR -FAR FIELD MEASUREMENT MODEL
The problem of estimating the direction to a transmitter, based on measuring the phase shift between antennas in a receiver array, is typically based on the measurement model [19] x(t) = A( , α)s(t) + n(t) (1) for a signal downconverted to baseband before sampling, where is the matrix of steering vectors for each signal, which uses the far-field assumption such that it depends only on the source direction parameterized by , α. This is a typical simplification that assumes that all array elements receive the incoming signal from the same direction as a plane wave with all elements sampled simultaneously. Defining the VOLUME 10, 2022 line-of-sight direction unit vector the steering vector, which is the predicted phase shifts for a given signal direction, takes the form where λ is the wavelength of the signal, and P a is the matrix of array antenna positions In our case, only a single transmitter is involved, and the model reduces to This model is not directly applicable for Bluetooth direction finding, because the signal sampled is not downconverted to baseband and the measurements are not simultaneous. This can be accounted for either by performing a measurement time correction on the data, as will be explained in Section III-B3, or by including an equivalent but opposite correction factor in the steering vector.
For an antenna array where elements are evenly spaced on a grid, calculating (4) with the full position matrix P a is not necessarily the most efficient, as the sine/cosine used in the calculation of the exponential may be slow depending on the computational platform architecture, and this is calculated for many direction pairs , α. An alternative method we propose is to calculate the phase step along each of the grid unit vectors, and then use complex multiplication to calculate the steering vector values. For the array in Fig. 1, we can calculate the phase difference between adjacent antennas with position differences p a x and p a y along the directions x a and y a , respectively, as Complex products can then be used to calculate the value for all elements where the index pair (n x , n y ) refers to the grid position of the antenna element along the unit vectors. To limit the number of operations, avoiding repeated multiplications, this can be performed stepwise as and so on for each direction, and combined as e.g.
Beamformer spatial pseudo-spectrum. α is 0 at the center, and 90 • at the edges, meaning that this covers the half-sphere in front of the array.

2) BEAMFORMER PSEUDO-SPECTRUM
The conventional beamformer [19], [20], also called the Bartlett beamformer, is a simple method for estimation of the direction from which a signal is received by an antenna array. This is essentially a spatial equivalent to the discrete Fourier transform, applicable to arbitrary array geometries. It is a spectral method where a spatial spectrum is to be searched to find the direction where the steering vector correlates best with the measurements. While algorithms such as Multiple Signal Classification (MUSIC) [21] and Estimation of Signal Parameters via Rotational Invariance Technique (ESPIRIT) [22] offer benefits such as increased resolution and better handling of multiple signals, beamformers have low computational complexity, without the need to compute the received signal correlation matrix or performing an eigenvalue decomposition. In our case, we only have a single transmitter and use a single-board embedded computer for real-time processing, and therefore use the beamformer approach.
For the measurement vector x and steering vector (4), the beamformer pseudo-spectrum is Fig. 5 shows an example of such spectrum, with dark red indicating the peak value. For Bluetooth CTE samples the measurement vector x typically contains multiple measurements from each array element. The steering vector must match the measurements, and if x contains more than one measurement from an element, the steering vector entry for that antenna must also be repeated. If the measurements have been corrected for CFO, such that they can be treated as being simultaneous, a general way to reduce the size of the steering and measurement vectors is to sum all measurements from the same antenna element, which correspond to the same complex steering vector value. The n k measurements x k,1 , . . . , x k,n k for antenna k is summed, and x = [x 1 , . . . , x m ] ∈ C m is used in (14). This yields the same estimation result as without the summation, as this only rearranges terms in the dot product of (14), Antennas with a higher number of measurements, and with higher measurement magnitude, will have a greater impact on the pseudo-spectrum. After summation a( , α) and x have only m elements without values for the same antennas repeated in a( , α), making the time spent on the spectrum peak search independent of the number of measurements for each antenna.
To find the estimated direction in a common navigation frame independent of the orientation of the array, we want to find the direction parameters n , α n . After finding the antenna frame direction pair , α maximizing the pseudospectrum, the line-of-sight vector l a ( , α) can be transformed to the local direction n , α n using where R n a depends on the array orientation. atan2(y, x) = Arg(x + jy) is the two-argument four-quadrant arctangent function standard in many programming languages, which is equivalent to tan −1 (y/x) for x > 0.

3) CFO ESTIMATION AND MEASUREMENT CORRECTION
To use (14) with the steering vector (4), we need to compensate for the CFO of the measurements. Since this can change over time it must be estimated, as shown in Section III-A. The unknown CFO needed to convert the measurements to baseband, allowing them to be treated as simultaneous, can be estimated using a similar approach as for direction, using the classical periodogram [18] adapted to nonuniform sampling and multiple antennas. A one-dimensional spectrum search is used to find the CFO estimate, if it is assumed constant throughout a CTE. If it is desirable to not assume a constant CFO, a linear drift term can additionally be estimated, making the spectrum two-dimensional, but this is not considered in the following. The Bluetooth specification [14] allows a maximum channel center frequency deviation of ±150kHz, corresponding to a clock error range of approximately ±60ppm, which comes in addition to the modulation deviation range mentioned in Section III. In the approach suggested here, each receiver antenna contributes individually to the frequency spectrum, and knowledge of the antenna geometry is not needed for the CFO estimation step. Any antenna element where multiple measurements are sampled during the CTE can be included, which helps to refine the frequency deviation estimate compared to only using the reference measurements.
Using only the reference period is the most straightforward and often suggested method for estimating the CFO [16], [23]. The low number of samples available from the reference period motivates the additional use of samples from the sample slots to improve the estimate. In [23] it is proposed to return to the same antenna every other sample, and include these repeated measurements in the frequency estimation. If the array has more than two elements, allocating half of the possible samples to one element obviously has the drawback of reducing the number of measurements available for direction estimation from the remaining elements. An analysis of several other approaches is found in [24], including both estimation of a global CFO for a whole CTE and local estimation where the frequency is considered to vary throughout the CTE. All methods considered there have frequency estimation ranges smaller than what the data allows with 1 µs spacing between reference samples, and smaller than the maximum possible CFO.
The frequency ''steering vector'' depends on the frequency and the time at which the included measurements were taken relative to each other, where t k is a column vector of the relative measurement times for antenna k, e.g. relative to the first measurement. A spectrum that can be used to estimate the CFO is then where x k is the subset of the measurement vector x with measurements from antenna k. An example spectrum using all 82 samples (8 reference samples and 74 from sample slots for 160 µs CTE with 1 µs slots) compared to the same method used with reference samples alone is shown in Fig. 6. Doppler shift due to movement of the transmitter or receiver along the direction between the two during the 160 µs CTE is absorbed into the frequency estimate. After finding the frequency with the maximum value of P freq (f ), this is used to correct either the IQ measurements themselves or the direction steering vector. The measurement correction is wheref is the estimated frequency deviation. After the correction, x = x 1 , . . . , x m can be used as if all measurements have the same time of validity. Example frequency estimates from a data set collected in the field, from extracting the peak VOLUME 10, 2022 FIGURE 6. Example periodogram frequency spectrum for a single CTE, from field data. The blue spectrum, with its scale on the left, uses all available measurements. The red, with its scale on the right, only uses the reference samples. The measurements correlate best with the steering vector for the frequency with the maximum spectrum value. Note the very different axis values; the estimate using only reference samples is much flatter than the spectrum using all available measurements.
frequencies in the spectrum for both cases shown in Fig. 6, is plotted in Fig. 7. The estimate using all measurements is also shown enlarged, and with filtering, in Fig. 8. The noise level is reduced considerably by the increased amount of measurements included.
The CFO estimate will be ambiguous with the frequency change required to change the phase of successive measurements by multiples of a full cycle. For a mix of different measurement intervals, sidelobe peaks will appear in the spectrum in addition to the ambiguity peaks, corresponding to each measurement interval. For reference measurements taken at the same antenna at 1µs intervals, this frequency ambiguity is (1µs) −1 = 1MHz. The ambiguities are not a problem for the CFO measurement correction as any of the equally strong peaks will produce the same corrected measurements.

4) PSEUDOSPECTRUM PEAK SEARCH STRATEGY
Once (14) is found from the measurements, a typical approach for finding the best fitting direction is to calculate the value P( , α) on a grid of α from 0 to 90 degrees and from 0 to 360 degrees with a chosen resolution. The direction estimate is then found as the direction with the highest value. The downside of this approach, which is well-known [25], is that the grid resolution must be high to get a small discretization error in the direction estimate, while at the same time keeping the resolution low enough to meet limits on computational resources, especially for real-time applications such as vehicle navigation. Multi-step searches using e.g. a coarse search followed by a fine search in a limited area around the coarse peak can make this more efficient, but it is still limited by the resolution of the final search.
Under the assumption that we are only receiving a signal from a single transmitter, and that any multipath either results in a combined peak, or a separate peak that is lower than the direct signal peak, we only need to find the highest value of the spectrum. This can be done efficiently using the following two-step process: 1) Run a coarse search. Knowing the geometry of the array, the expected power of any multipath peak compared to The few remaining small spikes in the estimates can be removed by outlier rejection, and noise can be lowered by low pass filtering, which is done in Fig. 8. the direct peak, and ideally the radiation pattern of the array elements, we can use the main lobe beamwidth to determine an appropriate search resolution that is as coarse as possible, while still ensuring that the highest spectrum value found is somewhere on the largest peak.
The resolution does not have to be uniform, as the main lobe size can vary for different search directions due to array geometry projection and the radiation pattern. Any directions that are considered impossible can be excluded from the search. As an example, the array power pattern of a square array, as shown in Fig. 1, but using omnidirectional elements, for a signal source in the boresight direction, was calculated using Matlab and is shown in Fig. 9. The essential point here is that a smaller array, which has a lower angular resolution, can use a coarser search than a larger array. 2) With the coarse direction found as the initial value, use a nonlinear programming (NLP) [26] solver to maximize P( , α). The spectrum is locally very smooth, especially for small arrays with a low angular resolution, and convergence to the maximum does not need many iterations. Nelder-Mead [26] gradient-free search can also be used, but for beamforming methods or MUSIC the spectrum functions are differentiable, meaning that the gradient can easily be computed, which should be taken advantage of. As the direction towards the transmitter is unlikely to change significantly between measurements (at e.g. 10 Hz measurement rate), it is not always necessary to run the coarse search for every new measurement set. Once a new measurement set is available the NLP solver can be used immediately, using the previous direction estimate or a prediction also based on other sensors, such as an inertial measurement unit (IMU) and an inertial navigation system, as the initial guess. Such sensors and systems are standard in many robotic applications.
If λ in the steering vector (4) is calculated based on the nominal carrier frequency, instead of the total frequency from the estimated CFO, the steering vector for all coarse search directions can be precomputed once and simply read from memory for each coarse search, which can speed up processing. The frequency deviation estimate is most important for measurement correction as explained in Section III-B3, while its use for wavelength calculation is of little importance as it only changes λ by single-digit parts per million.
If , α are used as the antenna-frame direction parameters, a uniform grid of these values will not be uniformly spread in direction, as illustrated in Fig. 10. Near the array boresight direction, neighboring azimuth values will be closer together. The result is that the search is inefficient, as some points are unnecessarily close. The azimuthal resolution close to boresight can be reduced as shown in Figs. 10b and 10c, or a different parameterization can be used as in Fig. 10d, to reduce processing time used for the coarse search. These are however simple examples that do not take into account knowledge of a specific antenna array.

IV. FIELD EXPERIMENTS
A. EQUIPMENT 1) GROUND ANTENNA EQUIPMENT Fig. 11 illustrates the hardware components connected to the antenna array. The array used is a Nordic Semiconductor experimental reference design, using 12 truncated corner right-handed circular polarization (RHCP) patch antennas in a square 15 × 15cm pattern, with 5cm antenna spacing. An nRF52833 board is directly attached to the back of the antenna using two header rows and a separate coaxial cable, where the header rows control the hardware signal switches in the array printed circuit board (PCB). The receiver board scans for advertising packets and outputs measurements only if the media access control address (MAC address) of the transmitter matches a pre-programmed value. Thus the system operates connectionless on the advertising channels with frequencies 2402 MHz, 2426 MHz and 2480 MHz. Due to the connectionless setup, Bluetooth was only used for direction estimation and not as a vehicle telemetry and command link, although this is a potential future extension. The board is configured to sample the CTE at a rate of 1 MHz, meaning once in every sample slot, for a total of 82 measurements.
The nRF52833 board is connected to a SentiBoard [27] universal asynchronous receiver-transmitter (UART) port. The SentiBoard forwards the data using Universal Serial Bus (USB) communication to a Beaglebone Black singleboard computer, where the data is parsed in DUNE [28]. A custom binary protocol is used for the output of the measurements from the nRF52833. This includes phase measurements as in-phase and quadrature components, the frequency of the channel used, and metadata about measurement timing and sampling order. The Beaglebone Black is either connected using an ethernet cable to the vehicle ground station, where a Ubiquiti Rocket M5 radio allows communication with the vehicle, or directly to a radio that is wirelessly connected to the ground station for remote array placement. The receiver data are transmitted to the payload computer of the vehicle where the direction estimation runs. A twostep coarse-fine spectrum search is used in these experiments, with coarse steps of 2.5 • and fine steps of 0.1 • . The CFO is estimated using 25 Hz search steps. A uBlox ZED-F9P GNSS receiver with a helix antenna is used as a Real-Time-Kinematic (RTK) GNSS base, transmitting measurements using the RTCM3 protocol to the UAV over the network, in order to evaluate the positioning performance of the Bluetooth system. The components are mounted to a custom bracket, pictured in Fig. 12, with slots at different angles for the array PCB relative to the base plate. This allows changing the pitch angle of the array while the bubble level on the bracket is used to level the assembly.    Fig. 13 shows a schematic with the relevant hardware components onboard the vehicle. A directional TrueRC Canada X-AIR 2.4GHz RHCP antenna is used. Fig. 14 shows the antenna mounted in the nose of a Skywalker X8 UAV pointing forwards. The nose of the UAV, in front of the antenna, is made of expanded polystyrene with a layer of canvas tape  for strength. Lab tests show no significant reduction in signal strength from this placement. The antenna is specified to have a gain of 8 dB, a −3 dB beamwidth of 75 • and performance equal to an omnidirectional antenna in a 120 • beam [29]. The antenna is connected to the nRF52833 board running transmitter firmware using a coaxial cable. The antenna connector on the board contains a signal switch that connects the transmitter to a linear polarization PCB trace antenna when an external antenna is not connected. A uBlox ZED-F9P GNSS receiver with a helix antenna is used on the UAV, receiving measurements from the antenna mounted on the array. The RTK GNSS setup yields very accurate and precise estimates of the UAV's position relative to the GNSS antenna attached to the array, with position errors on the centimeter level.

2) VEHICLE PAYLOAD
The transmitter broadcasts advertising packets with a 160 µs CTE at an average rate of 10 Hz (with a small random variation for conflict reduction with other Bluetooth devices). The SentiBoard outputs measurements from all connected sensors to the Odroid XU4 computer where they are both logged for later analysis and parsed for real-time use. The Odroid and SentiBoard are shown in Fig. 15, where they are mounted on top of the nRF52833.

B. PRACTICAL SAMPLING ASPECTS: SAMPLE RATE, SWITCHING INTERVALS AND CTE DURATION
The number of IQ measurements from a CTE is configurable in multiple ways. The CTE length can be configured, up to 160 µs, the duration of the switch and sample slots can be set to either 1 µs or 2 µs, and the sampling frequency can be set higher than 1 measurement in each slot. For the Nordic Semiconductor equipment used in the experiments in this paper, the maximum sample rate is 8 MHz, yielding 8 samples for each 1 µs sample slot, for a total maximum of 600 samples (or 1192 if sampling is also performed in the switch slots). Depending on the processing power available, the rate of CTE transmissions, and the algorithm used, this can be too much to process in real-time. Sampling can be done in both the switching and sampling periods and it can be configured in the external processing software whether any measurements in the last part of the switch slots are to be used for estimation if the switching transients settle early in the switch slots. A limitation of the communication interface between the nRF52833 board and the computer used for processing is the maximum UART baudrate of 1 Mb/s. With a start bit and stop bit, the data capacity is 800 kbit/s or an average of 80 kbit for each sample set at 10 Hz. With 32 bits being used for an IQ sample pair, 1 MHz sampling results in 82 IQ pairs for a total of 2624 bits, and 8 MHz sampling equivalently yields 19200 bits, or 38144 bits if switch slot samples are output as well, which is close to half the link capacity. If additional metadata is included, such as timestamps within the CTE or the number identifying the array element used for each measurement, or if the measurement rate is increased beyond 10 Hz, the required data rate could reach the link limit. Since the sample order and sample timestamp sequence within the CTE is always the same, they can be hard-coded in the processing software instead of output on the link, to reduce link usage and latency, although any changes in the receiver configuration would require changes in the external processing software.

1) MULTIROTOR
Initial testing was performed using a multirotor UAV, using the built-in linear polarization PCB trace antenna on the transmitter board, and the default 0 dBm transmit power. The array was placed flat on the ground pointing upwards, with axes aligned as best as possible to NED. The UAV was flown above the array at different heights up to 100 m, and at different horizontal distances. This is pictured in Fig. 16. Using height from RTK GNSS, Bluetooth direction measurements were transformed to North-East horizontal position for comparison with the RTK GNSS position estimates. Directly above the array, there was no significant loss of signal at 100 m height using the transmitter PCB trace antenna. The signal strength reduced with increasing horizontal distance, when the array received the signal at a more acute angle. Since the patch antenna elements in the array are somewhat directional, this is expected. The most interesting observation was that the direction estimation performance appeared significantly better for Bluetooth in the North direction as opposed to East, as clearly visible in Fig. 17 (the abbreviation BT is used for Bluetooth). This difference was independent of the UAV yaw angle, and when the array was yawed on the ground, the effect stayed the same in the antenna frame. Changing the array sampling order did not have any effect. Rotating the payload on the UAV by 90 degrees yielded the result in Fig. 18, which was much improved. Returning the payload to the     receiver onboard the UAV, which also uses 2.4GHz and sends telemetry data to the RC controller, in some way interferes with the Bluetooth system. For this reason, it was decided to try an external circularly polarized antenna, depicted in Fig. 19, which should be a better match for the circular polarization of the antenna elements on the array. This antenna is directional with a quite wide beam, and has a PCB ground plane that should shield the antenna from equipment above it on the UAV. Using the new transmitter antenna, Fig. 20 shows the result from a 20m diameter circle flown at 50m height, compared to RTK GNSS. The results for the three frequencies used are plotted separately in Fig. 21. Both directions appear similar in performance. It is visible that estimates from the three frequencies are slightly different, which could be related to a lack of frequency-dependent IQ origin calibration for the specific device. If the frequencies could be calibrated to remove the relative offsets, the performance looks very promising for multirotor precision landing without GNSS. Results for low elevation angles, α n < 40 • , were much worse than for high elevation angles, and the direction search in the spatial pseudo-spectrum was therefore limited to angles within 50 • of the array boresight, which was also carried over to the next fixed-wing test.

2) FIXED-WING
With the payload and RHCP antenna mounted in a Skywalker X8 fixed-wing UAV, tests were performed with the array  boresight direction close to horizontal, with a 10 • upwards pitch. Transmit power was increased to the maximum supported, 8 dBm, to provide the maximum range possible. Fig. 22 shows the array setup, which was leveled using a two-dimensional bubble level on the array mounting plate.
A flight was conducted where several maneuvers were tested, including straight lines towards the array at a constant height from multiple directions, crossing the boresight vector, changing both height and direction while approaching the array, and loiters at several distances, including one passing over the array. Figs. 23 and 24 show the RTK GNSS position of the UAV relative to the array, with color used to visualize the Bluetooth Received Signal Strength Indicator (RSSI) for all the Bluetooth measurements during this flight. The array was pointing South-West along the long straight line in Fig. 23, and Fig. 24 has a view from the side, perpendicular to the long straight line. The sensitivity of the receiver is −93 dBm according to the specification, and weaker signals will not be received. The maximum range in this test was approximately 700 m, which should be more than enough for the system to be usable in a UAV recovery scenario. An interesting observation is that at close range several spatial bands appear (see Fig. 25). Fig. 26 shows the azimuth and elevation angle estimates compared to RTK GNSS after calibrating the array heading by matching the azimuth angles from both systems.   The azimuth estimate from Bluetooth is very consistent, with no large systematic errors, although the array orientation assumed when transforming GNSS positions into the array frame can be slightly misaligned. The azimuth plots are shown enlarged in Fig. 27. The azimuthal standard deviation for the near-constant azimuth intervals 393-424s, 500-527s and 602-634s are 1.08 • , 1.00 • and 0.92 • respectively. The elevation on the other hand has very clear systematic and repeatable errors at elevation angles up to approximately 25 • , with larger errors for lower elevation angles. This was a trend for all maneuvers in this flight, with good azimuth performance and elevation angles struggling with systematic errors. The elevation angles are plotted separately for each channel in Fig. 28, showing systematic differences. To transform the Bluetooth angle estimates α, to horizontal and vertical antenna frame xy-positions, denoted p a xy,BT , the range along the array boresight vector from GNSS, −p a z,GNSS , was extracted from the NED GNSS position transformed to the array frame p a GNSS = R a n p n GNSS . Using p a xy,BT = −p a z,GNSS tan(α) we can see the position error the angles corresponds to in the middle part of Fig. 26. Fig. 29 shows corresponding plots for the four loiters with increasing distance shown on the left side of Fig. 23. An angular noise level will give rise to an increasing error in position as the range increases, but the angular noise level is decent at all distances where a signal was received. We expected direction estimation performance to become poor before losing the signal, so this exceeded our expectations.
The most likely source of the systematic errors in elevation is signal multipath. The signal received by the array is not only the signal propagating directly from the transmitter to the array but also a signal reflected from the flat ground surface, which causes signal interference and a signal that appears to originate from a different elevation angle.

V. MULTIPATH
If the signal broadcasted by the moving transmitter not only propagates directly to the receiver array but also reflects off surfaces and continues to reach the array, we have a multipath situation. For a long-range navigation scenario, the antenna array boresight may point close to horizontally, with the direction to the transmitter having a shallow elevation angle over the horizon. The signal transmitted can then be expected to reflect off the ground or water surface below, and the signal reaching the array will be the sum of a direct and a reflected signal, as illustrated in Fig. 30.
If the array is large enough, the angular resolution will allow the separation of the direct and reflected signals if the angle between them is sufficient. For a transmitter far away, the angle separating the direct and ground reflected signals is approximately two times the transmitter elevation measured from the antenna location, and separability thus requires a larger array for lower elevation angles. If the array is too small to satisfy the Rayleigh criterion for the transmitter direction and its reflection, a single combined peak appears in the spectrum. The Rayleigh criterion is satisfied if the angular separation between the peaks is at minimum as large as the separation of each peak and the first null of its response, known as the Rayleigh resolution limit [20]. For the array used here, Fig. 9 show this to be approximately 30 • if the first signal is in the boresight direction (90 • ), with the null at approximately 60 • along = 0. It is important to note that since the projection of the array element positions is smaller in other directions, the angular resolution of the array is reduced for directions away from the boresight. The interference of the two signals will result in an error in the direction of the pseudo-spectrum peak compared to the pseudo-spectrum of only a direct signal.
Signal reflection behaves differently for different frequencies, surface materials, and reflection angles. For a circularly polarized signal, the reflected signal is linearly polarized if the angle of incidence equals Brewster's angle, where θ B is measured from the surface normal. n surface and n air are the refractive indices of the surface and air, respectively. A RHCP signal reflected off a surface at normal (i.e. 0 • ) incidence will have left-handed circular polarization (LHCP). For incidence at the Brewster angle θ B , the reflected signal is linearly polarized in the horizontal direction [30]. For angles between the Brewster angle and the normal, the polarization is left-handed elliptical. For signals incident at large angles of incidence, or small grazing angles, the reflection will remain right-handed polarized, but also elliptical.
Since an RHCP receiver antenna rejects LHCP signals, incidence at near-normal angles results in less multipath issues than small grazing angles, which is the case for the experimental results in Fig. 26.
Since the refractive index of the surface material is higher than that of air, the phase of the reflected signal changes by a half cycle. For a circularly polarized signal, this occurs for all incidence angles [30].

A. MULTIPATH SIMULATION
A simulation study was conducted to determine how a ground reflection would manifest in direction estimates. Assumptions are a completely flat ground, and a half cycle phase shift in the reflection due to refractive index. The distance traveled by the direct signal path between array element i and the transmitter antenna is where p n a,i is the North-East-Down (NED) position of the i-th receiver array antenna element, and p n t is the position of the transmitter antenna. The complex phase measurement of receiver antenna i is calculated from the distance as Note that the measurements for each antenna, unlike the real measurements, are simulated as simultaneous. For a completely flat ground surface, assuming that the angles of incidence and reflection are equal as illustrated in Fig. 30, the coordinate of the reflection point p r,i for the i-th array element is   where subscript NE denotes the two-dimensional North-East component of the vector and subscript D denotes the scalar Down (vertical) component. The path distance traveled by the reflected signal is and the phase measurement for the reflected signal alone is A coefficient c r , typically between 0 and 1, is used for attenuation of the reflected signal. It should be noted that the phase VOLUME 10, 2022  is inverted by the addition of πj in the exponent. The phase measured at the receiver antenna is considered to be the sum of the direct and reflected signals, These simulated measurements can then be used in the direction estimation algorithm in the same way as real measurements but without any need for frequency estimation or correction for the difference in measurement time. To visualize the error in estimated elevation angles for different positions in front of the array, we estimate direction for positions on a grid of heights and horizontal distances from the array. Simulated measurements are created for a transmitter at each position, and the pseudo-spectrum is searched using an NLP solver initialized using the known true direction. Simulation results are shown in Figs. 32b to 35b compared with experimental data explained in Section V-B.

B. EXPERIMENTAL MULTIPATH TESTING
To validate that the effect of multipath on the elevation angle estimate behaves as described in Section V-A, a field experiment was performed using a DJI S1000 multirotor UAV carrying the Bluetooth transmitter board and antenna. The test was performed on a grass runway covered with snow, shown in Fig. 31. The antenna array was set up on three different heights, with the array center at 1.35 m, 0.7 m and 0.15 m above ground, and upwards pitch angle of 10 • . The mounting bracket had a 10 • angle built-in and was leveled using a two-dimensional bubble level. Finally, a test was conducted at roughly 0.17 m height with a pitch angle of approximately 45 • without leveling (later estimated as 48.5 • degrees from the data). The runway was quite flat, see Fig. 31, but the slope was not perfectly known. RTK GNSS was used for accurate and precise relative positioning, for comparison with Bluetooth estimates, as in the previous experiments. The UAV was flow on a flight path along the runway in front of the array, with approximately a constant azimuth angle n , distances out to 180m and a height up to 100m. The flight path consisted of straight line segments, as shown using a side-on view in Fig. 32 where the array is located at the origin in the bottom right corner. The error in the Bluetooth elevation angle estimates were computed by subtracting the elevation calculated from RTK GNSS position measurements, which were used as ground truth values. Figs. 32a to 35a shows the computed Bluetooth elevation angle errors from the field experiments, using only the 2480 MHz channel, using color along the flight trajectory. Red color indicates that the Bluetooth elevation estimate is higher than the RTK GNSS elevation, and the opposite for blue color. Figs. 32b to 35b shows the predicted result from simulated measurements for the same array placement scenarios as in the experiments. For a position grid of distances and heights, elevation angles were estimated using simulated measurements with both c r = 0.25 and c r = 0 (multipath-free). By subtracting the multipath-free estimates from the estimates with multipath, the simulated error was found. The dotted lines are drawn with the same angle in both figures, centered in positive error rays from the experimental data, for easier comparison of the patterns.
The simulation does not take polarization of the reflected signal, or antenna properties such as radiation pattern, into account. The polarization mostly causes disagreement at higher elevation values, where the error is reduced since the reflection is then left-hand polarized, which is not what the array is designed to receive. The effect of radiation pattern has the most effect in the scenario in Fig. 35 with the array pointing upwards in the direction of the UAV when flying at high elevation angles. The reflection is then received by the array from a direction almost in the array PCB plane, where the gain is much lower for the directional patch antennas.
Comparing the experimental and simulated results, where the experimental results are considered samples of an underlying error pattern covering the area shown in the simulated results, it is clear that the observed pattern is similar to the simulation, especially for low elevation angles. If the offset between the experimental and simulated results was the same for all three heights, a reasonable source of the offset would be runway slope or a bias in the array leveling. Since the array was re-leveled for each of the heights over the ground, differences in the residual leveling error would lead to differences in the resulting angle offsets. The location on the ground in front of the array where signals reflect depends on the array height above ground, being closer to the array for low array heights. The simulation assumes a perfectly flat ground, which is not the case in reality, being a significant error source.

C. FACTORS AFFECTING ELEVATION ERROR
Based on simulating different scenarios, the effect of a number of parameters was observed, which are shown in Fig. 38 compared to the reference case Fig. 38a: Fig. 38b) A stronger reflection increases the error amplitude. This depends on ground material including water content, e.g. wet grass reflects more than dry grass. Seawater is highly reflective for radio frequencies.
2) HEIGHT ABOVE SURFACE (Fig. 38c) Increasing the height of the array over the surface increases the spatial frequency of the elevation error oscillation with respect to the true elevation angle. To explain why these oscillations occur, we combine (14) with the multipath signal sum (30), resulting in For elevation angles where the complex value a( , α) x reflected aligns with a( , α) x direct in the same or opposite direction, the addition of the reflected signal does not influence the estimated elevation angle, and the error is zero. When the components do not align, the elevation angle maximizing the spectrum changes. Note that the wide ''bands'' of reduced error remain in the same directions, and that these depend on the amplitude of the error, not the phase. Minima in the amplitude occurs when the argument (phase angle) of the values of the complex vector a( , α) x reflected spread out, such that their sum a( , α) x reflected has a small magnitude. On the other hand, when the values align such that a( , α) x reflected has a maximum magnitude, the error amplitude is maximized.
3) ANTENNA PITCH ANGLE (Fig. 38d) By pitching the array up there is essentially a reduction in the vertical array size seen from the front through projection. In the simulated case this moves zero-amplitude rays to higher angles without changing the phase of the oscillations.
A second effect that is not included in the simulation is that by pitching the array, the radiation pattern maximum, which FIGURE 34. Multirotor validation test using RTK GNSS, array center 0.15 m over ground and 10 • upwards array pitch. In this case, the upper dashed line was placed in a positive error peak for the experimental data, but it appears that the minimum-amplitude direction between the dashed lines for simulation occurs at a higher elevation value for the experimental data, resulting in a very narrow positive peak at the edge of the red beam visible in the simulation result. The increased negative error between the dashed lines in the experimental data also supports this.

FIGURE 35.
Multirotor validation test using RTK GNSS, array center approximately 0.17 m over ground and 48.5 • upwards array pitch. Compared to 10 • tilt, the 25 • zero amplitude lobe is moved up to approximately 40 • , and the oscillation spatial frequency is increased a little for low angles due to the effectively smaller array from this direction. points in the boresight direction, points upwards. This reduces the gain in downwards directions from which the reflected signals are received. Additionally, the directions blocked by the array ground plane change. Consider a pitch angle where the direct signal has 0 dB gain while the reflection has −3dB, and compare this to a pitch angle where the direct signal has gain −3dB while the reflection has −10dB. The effect of multipath will be reduced, at the cost of reduced range. This can go as far as having the reflected signal arrive at the array from behind.  increases their frequency. This reduces multipath error, but also requires an increased number of elements to avoid ambiguity errors or too high sidelobe power. The orientation of the array also matters if the maximum dimension is different in different in-plane directions, e.g. a square antenna can be mounted with a diagonal vertically, giving an increased vertical size, as illustrated in Fig. 36.
For a signal source at an infinite distance (far field) in front and higher than the array, the reflected signal will reach the array at an angle from below equivalent to the angle from above, illustrated in Fig. 37.  (Fig. 38g) For Bluetooth, the carrier frequencies of the different channels give a minor effect on the error, which is more pronounced at higher elevation angles. If higher frequencies could be used, however, such as 5.8 GHz, multipath would be reduced since an array of the same physical size would be larger relative to the signal wavelength. This would however require more elements to avoid ambiguity errors.

D. CALIBRATION
Canceling the multipath elevation error using a calibration for a specific scenario could be beneficial for the practical use of Bluetooth for navigation. A simple elevation angle calibration is however not globally possible above a reflection threshold that depends on the situation, as the true and estimated elevation angles do not always have a one-to-one relationship: More than one geometric elevation angle can result in the same measurement due to the multipath error. Fig. 39 shows simulated examples of estimated elevation angle plotted as a function of the true elevation angle for a transmitter at 500 m distance from the array, with the different colors representing different reflection coefficients. This is created using the method described in Section V-A, but with positions at a constant Euclidean distance with a range of elevation angles. One way the calibration problem could be handled is to track the direction over time, but this does not solve the problem if the initial direction is unknown. Fig. 40 shows a similar plot for field experiment data.

VI. CONCLUSION AND FUTURE WORK
The use of Bluetooth direction finding for outdoor navigation has been shown to be feasible at distances up to 700m, close to the receiver sensitivity of -93dBm, without significant loss of angular precision, which exceeded our expectations. Azimuth angle estimates were accurate, with standard deviation on the 1 • level, while the elevation estimate suffered from systematic errors due to multipath. The effect of multipath was investigated using simulation and field experiments, showing quite consistent results. The range and performance of the system look promising for fixed-wing UAV recovery and ship docking in GNSS-denied environments, which will be a focus of future work. The proposed peak-search method should be implemented for online use in a payload, benchmarking its computational requirements compared with other methods. Methods to reduce the effect of multipath should be considered.