# Low Cost Antenna Array Based Drone Tracking Device for Outdoor Environments

Authors: Marcos T. de Oliveira; Ricardo K. Miranda; João Paulo C. L. da Costa; André L. F. de Almeida; Rafael T. de Sousa
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

Applications of direction of arrival (DoA) techniques have dramatically increased in various areas ranging from the traditional wireless communication systems and rescue operations to GNSS systems and drone tracking. Particularly, police forces and security companies have drawn their attention to drone tracking devices, in order to provide the safeness of citizens and of clients, respectively. In this paper, we propose a low cost antenna array based drone tracking device for outdoor environments. The proposed solution is divided into hardware and software parts. The hardware part of the proposed device is based on off-the-shelf components such as an omnidirectional antenna array, a 4-channel software defined radio (SDR) platform with carrier frequency ranging from 70 MHz to 6 GHz, a FPGA motherboard, and a laptop. The software part includes algorithms for calibration, model order selection (MOS), and DoA estimation, including specific preprocessing steps and a tensor-based estimator to increase the DoA accuracy. We evaluate the performance of our proposed low cost solution in outdoor scenarios. According to our measurement campaigns, we show that when the array is in the front fire position, i.e., with a DoA ranging from-60° to 60°, the maximum and the average DoA errors are 6° and 1,9°, respectively.


---

## Body

## 1. Introduction

Applications of direction of arrival (DoA) techniques have dramatically increased in various areas ranging from the traditional wireless communication systems [1, 2] and rescue operations [3] to GNSS systems [4–7] and drone tracking in public and private events. In the last years Unmanned Aerial Vehicles (UAVs) have been a major concern of airspace control bodies and military due to possible terrorist attacks and illegal activities. In 2015, there were more than nine hundred incidents involving drones and aircrafts in the United States [8], whereas, in April 2016, a UAV reached an aircraft landing at the Heathrow airport in London [9]. In 2016 in Dubai, four drones invaded the airport interrupting the landings and take-offs, causing an estimated loss of one million dollars [10]. In October 2017 in Canada, the first reported collision of a drone and a commercial airplane has occurred [11]. Recently police forces and security companies have drawn their attention to drone tracking devices in order to provide the safeness of citizens and of clients, respectively. In this sense, the development of low cost devices for drone tracking is fundamental to fit such demands.In general, the DoA estimation techniques can be broadly classified into conventional beamforming techniques, maximum likelihood techniques, and subspace-based techniques [12–14]. In [15] the authors proposed to estimate the DoA of a signal impinging the Electronically Steerable Parasitic Array Radiator (ESPAR) antenna with twelve parasitic elements by using a support vector machine (SVM) technique. In the anechoic chamber the result of the experiments reaches 0.67° estimation error. No hardware details were provided.To overcome the effects of multipath propagation on the performance of DoA estimation, the authors of [16] proposed a frequency domain multipath resolution subspace-based approach, which makes RSS-based DoA estimation applicable in multipath scenarios for small-size and low-power sensor networks. This approach was verified with Monte Carlo simulations with high SNR.Aiming at the growth of connected cars systems, the authors in [17] developed a 4×4 MIMO antenna system and proposed the DoA function for a circular phased array antenna. Since the focus of this paper is on the methodology and basic characteristics of DoA function, no measurement or hardware information was provided.In [18], DoA estimation using an ESPAR with 12 parasitic elements and one active monopole is carried out for wireless sensor network (WSN) applications. The authors calibrated the ESPAR array using an anechoic chamber. Since the focus of [18] is on the calibration, no outdoor or indoor measurement campaigns were performed by the authors.An improvement of ESPAR antennas is proposed in [19]. The authors developed a Multiple Beam Parasitic Array Radiator (MBPAR) antenna that can realize six beams at the same time without the use of diodes, which increases the communication capacity. To validate the proposed design, a prototype was fabricated at 2.45 GHz. The antenna has the efficiency from 94.2 % to 95.7 % over the 2.4G-WLAN bands. No DoA measurements were provided.In [20], a square shaped 16 element antenna array is connected to switches so that a four-channel SDR can select four antennas at each side of the square, allowing a 360° DoA estimation in outdoor environments. Each side of the square performs a ±45° azimuth estimation. According to the authors, three Yagi antennas were used as sources at specific points, and a maximum DoA error of 5° is achieved. No information is provided about the real distance between the sources and the receive array.In [21], a four-element quasi-Yagi antenna array system is applied for DoA estimation using the MUSIC algorithm, whereas the Minimum Description Length (MDL) criterion is used to estimate the number of dominant multipath components. Only two measurements were performed for two specific positions, showing an error of 1°. However, no information is provided about the experimental scenario.In [22], several DoA estimation techniques are compared considering a horizontal uniform linear array (ULA) with 12 elements inside an anechoic chamber. The measurements were conducted varying the DoA from -20° to 20° in steps of 4°. The DoA estimation errors were smaller than 2°. The Min-Norm approach [12] outperformed MUSIC [23] although it has a higher standard deviation.Finally, in [24], the authors developed system using five-port reflectometers that allow simultaneously measuring the DoA and Time of Arrival (ToA) of coherent and incoherent signals, connected to seven quasi-Yagi antennas, with one reflectometer for each antenna. The MUSIC algorithm is applied for the DoA estimation, providing an error of 2° for one source and 0,5° for two sources. The measurements were performed in a nonreflective environment.In order to detect the presence of drones and to track them, there is a variety of mechanical, optical, or antenna array based solutions in the market. For instance, the mechanical solution in [25] detects a drone within 3 km for targets up to 55 cm in diameter and classifies the model of the drone within 1.1 km. The position accuracy (azimuth) in [25] is 1°. In [26], a rechargeable portable drone tracking device can detect and indicate the direction of a drone in a 360° plane even with weak line of sight (LoS) component. The device in [26] allows the communication with other devices by using an Application Programming Interface (API) framework. No technical information and patent about the principles behind the device in [26] and its DoA accuracy were provided. In [27], an antenna array based system is shown to detect with a 1 km range and with 1° accuracy or with a 7 km range and with 3° accuracy.In [28], a mobile application (app) is proposed for drone detection. According to the developers, the app has an average range of 106 meters. The system allows the detection of almost 95 % of all types of drones. However, the solution in [28] does not indicate the position or the direction of the drone.In terms of drone tracking, there are recent works on received signal strength (RSS)-based DoA estimation. For example, [29, 30] propose to estimate the DoA using arrays of Yagi-Uda directional antennas for the localization of drones exploiting by their incoming NTSC signal in a measurement campaign. The work [29] proposes a complete hardware and software framework using arrays of directional antennas and formulates a novel DoA estimation correction procedure. In [30], a novel DoA estimation algorithm for the localization of drones is validated by using an AD-FMCOMMS5-EBZ software defined radio (SDR). Finally, [31] implemented tests to detect and locate UAVs at 900MHz. The authors used MUSIC combined with spatial smoothing and MDL. The work [30] has a similar objective to our proposal; however the authors did not concern in to present the accuracy of the system and in use of low cost equipment.In this paper, we propose a low cost antenna array based drone tracking device for outdoor environments. To the best of our knowledge, there are no state-of-the-art low cost off-the-shelf antenna array based devices applied to drone tracking. The problem of drone tracking is challenging due to the several possible modulation schemes for the data transmission, multipath propagation, and the possible long operational distances. The proposed framework can also exploit tensor-based techniques, such as the Parallel Factor Analysis (PARAFAC). In contrast to the subspace-based methods, the tensor-based approach shows to be robust in real scenarios.The remainder of this paper is divided as follows. In Section2, the data model is presented. Next, in Section 3, we propose a low cost antenna array based drone tracking device for outdoor environments, including a complete description of the hardware and software, and the steps involved for assembling, calibration, and signal processing. In Section 4, we validate our proposed solution by means of measurement campaigns in an outdoor scenario. In Section 5, conclusions are drawn.


## 2. Data Model

We assumed be far-field sources transmitting narrow-band signals. These planar wavefront signals impinge over a receive antenna array with M omnidirectional elements that are uniformly and linearly disposed. The space Δ between two adjacent antennas is equal to λ/2, where λ is the wavelength of the carrier signal. The received signals at the antenna array can be written in a matrix form as follows:(1)X=AS+N∈CM×N,whereA∈CM×d is the steering matrix and its i-th steering vector is given by(2)aμi=1ejμie2jμi⋯ejM-1μiT∈CM×1,where μi is the spatial frequency that can be mapped into the direction of arrival of the i-th source, θi, by the following expression: μi=2πΔsin⁡θi/λ.S∈Cd×N is the symbol matrix withN being the number of snapshots.N∈CM×N stands for the noise matrix whose elements are assumed to be Complex-Valued Circularly Symmetric Gaussian and identically and independently distributed (i.i.d.) random variables.Given (1) and assuming that the noise and the signal are uncorrelated, the covariance matrix can be computed by(3)RXX=ExxH=ARSSAH+RNN,where x is one column vector fromX, (·)H is the Hermitian operator, and E· is the expected value operator. In practice, the sample covariance matrix is calculated as follows:(4)R^XX=XXHN∈CM×M.The DoA techniques used along this paper exploit the sample covariance matrix in (4). As shown in Section 3, the matrixX is preprocessed before we compute the sample covariance matrix R^XX.The goal of our proposed drone tracking device is to estimate the direction of arrival (DoA)θ1 of the line of sight (LoS) component from a drone in an outdoor scenario. We assume that there is no obstruction of the LoS component. Therefore, the LoS component is assumed to have the greatest power in comparison with the non-LoS components. Mathematically, we can express it as(5)aμ1sμ1F>aμisμiF,for i=2,…,d. The operator ·F stands for the Frobenius norm.


## 3. Proposed Low Cost Antenna Array Based Drone Tracking Device

In this section, we detail the proposed low cost antenna array based drone tracking device. In Section3.1, we describe the steps for the hardware calibration. The calibration ensures phase alignment for all the four channels of the SDR, allowing the DoA estimation. In Section 3.2, we present the assembling of the hardware components of the proposed drone tracking device. In Section 3.3, we propose a signal processing framework for DoA estimation.


### 3.1. Hardware Assembling for the Calibration

In order to perform the measurements, the SDR should be calibrated, such that all the receive channels become in phase. The phase imbalance may be caused by different time initialization of the oscillators and by hardware imperfections. The hardware vendor provides a software [32] for clock calibration of the local oscillator. However, this software does not perform phase calibration.Therefore, in order to perform the phase calibration, the hardware components are first assembled according to Figure1. Note that the SDR transmits the signal from one channel and receives it in four channels that should be calibrated.Figure 1
Assembled components for the hardware calibration of four receive channels. The components are a microprocessor (uProc), a FPGA motherboard, a SDR, a power divider, and cables.As shown in Figure1, the SDR is a 4×4 MIMO platform named ADFMCOMMS5 [33], with two AD9361 [34] Integrated Circuits (ICs) that contain 2 transmitters and two receivers each, ranging from 70 MHz to 6.0 GHz, and have a channel bandwidth ranging from 200 kHz to 56 MHz. The platform is connected to a microprocessor and a FPGA motherboard [35] that configures the SDR and transmits the SDR data to the PC. As shown in Figure 1, the cables for calibration should have the same length. Moreover, a power division component is included in order to lead the signal to the four receive channels at the same time and to reduce the power of the transmitted signal to avoid damaging.To compensate the phase errors previously explained, the first step is to extract the phase of the elements of the matrixX. The phase ϕ(m,i) is defined as follows:(6)ϕm,i=∠xm,i=arctan⁡Im⁡xm,iRe⁡xm,i,where x(m,i) is the element in position m,i of the measured matrix X. The operators ∠, Im⁡· and Re⁡· stand for the phase operator, the imaginary part of a complex number and the real part of a complex number, respectively.Since each channel is related to each line of the matrixX, in order to compute the phase shift between two channels, we compute the phase difference of two consecutive antennas(7)ωm,:=Φm,:-Φr,:∈C1×N,where Φ stands for the matrix containing the calculated phases by (6), r indicates the reference channel, and m varies from 1 to M. This reference channel can be randomly selected from 1 to M and is the input of the SDR that is used as a reference to compensate the phase imbalance from the other inputs. Since the vector ω(m,:) is the m-th row of matrix Ω, in case m=i, the i-th row of Ω is filled with zeros. Finally, since the phase difference may slightly vary for different samples in the same row of Ω due to the thermal noise, we compute the arithmetic mean of the elements of each row of Ω, obtaining the vector ω¯∈CM×1 and its m-th element is given by(8)ω¯m=1N∑n=1NΩm,n.Hence, in order to compensate the phase shift between two different channels, we compute the vectorc∈CM×1. The m-th element is given by(9)cm=e-jω¯m.Note that the compensation vectorc is computed only once for the system initialization. The calibrated outputs of the antenna array are given by the following expression:(10)Xc=diag⁡cX,where the operator diag· transforms its argument vector into the main diagonal of a diagonal matrix.


### 3.2. Hardware Assembling for the Drone Tracking Measurement Campaign

After the hardware has been calibrated, the next step is to assembly it in order to perform the measurements.The four-element omnidirectional antenna array is connected to the calibrated hardware composed of the FPGA motherboard and SDR daughterboard according to Figure2. Each antenna is dual band (from 2400 MHz to 2483.5 MHz and from 4900 MHz to 5875 MHz ) [36] and has linear polarization with 3.7 dBi of gain. The space Δ between two consecutive antennas is equal to 5.99 cm. The operational frequency f=2.5GHz is the maximum frequency that avoids aliasing.Figure 2
Proposed low cost antenna array based drone tracking device with 4 element ULA and a 4× 4 MIMO SDR.


### 3.3. Framework for DoA Estimation

Here, we first propose a sample selection approach for DoA estimation by automatic phase deviation detection. Then, we formulate a DoA estimation framework exploiting preprocessing techniques and model order selection schemes.Figure3 depicts the flowchart of the proposed signal processing solution for DoA estimation.Figure 3
Flowchart of the proposed solution for DoA estimation.As shown in Box 2 of the Figure3, the phase deviation correction proposed in Section 3.1 returns a matrix Xc that is used in the sample selection scheme in Section 3.3.1.


#### 3.3.1. Sample Selection for DoA Estimation by Automatic Phase Deviation Detection

As exemplifiedXc in Figure 4, empirically we observed that the hardware causes phase deviations on the samples in random time instants. Therefore, we propose an approach to select the samples with phase deviations for the DoA estimation.Figure 4
Phase of the received data in each antenna. An example of phase deviation region is highlighted in the zoom area.Note that the phase compensation proposed in Section3.1 has been applied to the samples, whose matrix Φ containing the phases are depicted in Figure 4. Furthermore, note that there are significant deviations that can degrade the DoA estimation process. The main objective here is to remove these phase deviations.As shown in Figure5, such ripples can be better visualized by computing the phase difference in the time dimension according to the following expression:(11)γm,i=ϕm,i+1-ϕm,i2,where γ(m,i) is the value containing the quadratic difference of two consecutive time samples i and i+1 of the m-th channel. The ϕ(m,i), from (6), is the element in position (m,i) of the matrix Φ.Figure 5
Difference between consecutive samples from Figure4.Figure5 draws the Γ that contains the result of (11). By detecting the peaks, we can identify which samples should be removed. For this task, we apply the approach proposed in [37], which returns the green curve with the value of the threshold. Therefore, the samples whose phase differences are greater than the threshold are removed. The result after the samples removed is presented in the following equation:(12)Xcs=Xc:,1:N1∣Xc:,N2:N3∣⋯∣Xc:,NT-1:NT,where Xcs is the matrix with the selected samples. Note that NT=N and the values of Nt, for t=1,…,T, are found by comparing the phase difference values with the threshold in Figure 5.


#### 3.3.2. DoA Estimation Framework

According to Figure3, the matrix Xcs given in (12) is used to improve the DoA estimation with preprocessing schemes. There are several DoA estimation schemes in the literature, such as beamforming approaches and subspace-based approaches. Examples of beamforming are Delay and Sum [38] and Capon [39], whereas examples of classical subspace-based approaches are MUSIC [23] and ESPRIT [40].The DoA estimation schemes assume that the model orderd is known. In practice, model order selection techniques should be applied to estimate the model order d, as depicted in Figure 3.In the flowchart of Figure3, we adopted the Exponential Fitting Test (EFT) [43, 44] as the model order selection scheme. The EFT has the deflation property that allows us to find suitable thresholds as a function of the Probability of False Alarm (Pfa). By exploiting the deflation property and by finding suitable thresholds, the EFT has been the only scheme in the literature to estimate d=1 in the presence of a strong LoS signal and d=0 in the only noise (no signal) measurements. We have compared several schemes in the literature such as Akaike Information Criterion (AIC) [45], Efficient Detection Criterion (EDC) [46], Minimum Description Length (MDL) [47], Stein’s unbiased risk estimate (SURE) [48], RADOI [49], ESTimation ERror (ESTER) [50], and Subspace-based Automatic Model Order Selection (SAMOS) [51]. The M-EFT [43, 44, 52] has also been suitable, but an even smaller Pfa was required to find the thresholds. The computation of the thresholds of the EFT requires an extremely low Pfa. The complexity of such a computation is prohibitive. Therefore, we propose in Appendix A an extrapolation algorithm to compute such thresholds. Note that our proposed extrapolation algorithm has been applied in [53–56], although no details are provided. The reason for extremely low Pfa may be related to the noise behaviour as shown in Appendix B. Note that the Ilmenau Package for Model Order Selection (IPM) [57] with MATLAB and Java implementation of the model order selection schemes can be found at the LASP homepage (https://lasp.unb.br/index.php/publications/softwares/).In order to further improve the accuracy of DoA estimation schemes, preprocessing schemes can be applied beforehand. We consider in this work the Vandermonde Invariation Technique (VIT) [58], Spatial smoothing (SPS) [59], and Forward Backward Averaging (FBA) [60, 61] as a preprocessing schemes. As depicted in Figure 3, after the preprocessing step, a matrixZ is returned and used by the DoA methods summarized in Table 1.Table 1
Selected state-of-the-art DoA estimation schemes.
Delay And Sum [38] P D S θ = w θ H R x x w θ w θ H w θ . ( 13 )   Capon [39] P C A P ( θ ) = 1 w θ H R x x - 1 w θ . ( 14 )    MUSIC [23] P M U S I C ( θ ) = 1 w H θ V n V n H w θ, ( 15 )   ESPRIT [40] Ψ = J 1 U s + J 2 U s ( 16 ) Ψ = E V D E Φ E H ,   w i t h   Φ = d i a g  [ ϕ 1 , … , ϕ d ] ( 17 ) μ i = ∠ ( ϕ i ) ,   i = 1 , … , d ( 18 )In Table1, the vectorw(θ) in (13), (14), and (15) vary according to the candidate values of θ. The value of θ that maximizes the expression in (13), (14), and (15) is the θ^1, since, in Section 2, the data model assumes that the LoS component faces no obstacles. Therefore, the component corresponding to the greatest power should be the same component with DoA equal to θ^1. Us∈CM×d is signal subspace, which is equal to the d eigenvectors corresponding to the d greatest eigenvalues, whereas Vn∈CM×M-d is a basis for the noise subspace, composed by the M-d eigenvectors associated with to the M-d smallest eigenvalues. In (17) and (18), Φ is the diagonal matrix that has the eigenvalues of Ψ. We compute all the spatial frequencies and the one whose component has the greatest power is the μ^1 that can be mapped to θ^1.


### 3.4. Tensor-Based DoA Estimation

In this subsection, a tensor factorization, namely, the PARAFAC decomposition, is applied. The PARAFAC decomposition generates three factor matrices from a received tensorX, whose structure is detailed in this subsection. One factor matrix corresponds to the estimate of the steering matrix A containing DoA information. In this subsection A is extracted from tensor X to estimate the DoA of the impinging signal.We first consider an unchanging sequence ofN symbols transmitted periodically. Such symbols can be found in a header or footer or even in the payload of a message. Alternatively, repeating sequences of symbols can be extracted from time periods when no data is being transmitted, but the carrier of the transmitter is active. An example is given as follows. Since the oscillators at the transmitter and at the receiver are never exactly the same, a small frequency deviation or constant phase change is observed at the receiver. At an MSK receiver, if the deviation is positive, a sequence of ones can be extracted and if the deviation is negative a sequence of zeros is observed.The symbols corresponding to the sequence from thei-th source is represented by the vector ci∈CN×1. Accordingly, we can build a received signal matrix19Xip=γipaiμiciT,∈CM×Nwhere p is the period corresponding to the sequence transmission. Generalizing for d signals we find20Xp=ADγpCT,∈CM×N,where Dγ(p)=diag(γ1(p)γ2(p)…γd(p))T and C=[c1c2…cd].ForP transmitted sequences, the symbols can be concatenated along the third dimension to form the received signal tensor21X=X1⊔3X2⊔3…⊔3XP∈CM×N×P,where ⊔3 represents concatenation along the third dimension. Since the slices of X can be written as 20, X has a PARAFAC structure and can be decomposed into three factor matrices A, C and Γ, where Γ contains the diagonals of Dγ(p) along its rows. To factorize X, we first rewrite it in three different matrix representations or unfoldings22X1=AC⋄ΓT,23X2=CΓ⋄AT,24X3=ΓA⋄CT,where ⋄ is the Khatri-Rao (column-wise Kroenecker) product.It is known that minimizing22, 23, and 24 in the least squares sense [62] leads to the following solutions:25A^=X1C⋄ΓT+,26Γ^=X2Γ⋄AT+,27C^=X3A⋄CT+.We consider the well-known Alternating Least Squares (ALS) algorithm to solve25, 26, and 27 in an iterative way. Since C is known, the estimation step 27 is skipped, and the ALS algorithm alternates between the estimations of A and Γ in a two-step approach [63].Once convergence is achieved, we useA^ to extract the DoA of the i-th source as follows: 28θ^i=argmaxθ⁡aiHθa^i


## 3.1. Hardware Assembling for the Calibration

In order to perform the measurements, the SDR should be calibrated, such that all the receive channels become in phase. The phase imbalance may be caused by different time initialization of the oscillators and by hardware imperfections. The hardware vendor provides a software [32] for clock calibration of the local oscillator. However, this software does not perform phase calibration.Therefore, in order to perform the phase calibration, the hardware components are first assembled according to Figure1. Note that the SDR transmits the signal from one channel and receives it in four channels that should be calibrated.Figure 1
Assembled components for the hardware calibration of four receive channels. The components are a microprocessor (uProc), a FPGA motherboard, a SDR, a power divider, and cables.As shown in Figure1, the SDR is a 4×4 MIMO platform named ADFMCOMMS5 [33], with two AD9361 [34] Integrated Circuits (ICs) that contain 2 transmitters and two receivers each, ranging from 70 MHz to 6.0 GHz, and have a channel bandwidth ranging from 200 kHz to 56 MHz. The platform is connected to a microprocessor and a FPGA motherboard [35] that configures the SDR and transmits the SDR data to the PC. As shown in Figure 1, the cables for calibration should have the same length. Moreover, a power division component is included in order to lead the signal to the four receive channels at the same time and to reduce the power of the transmitted signal to avoid damaging.To compensate the phase errors previously explained, the first step is to extract the phase of the elements of the matrixX. The phase ϕ(m,i) is defined as follows:(6)ϕm,i=∠xm,i=arctan⁡Im⁡xm,iRe⁡xm,i,where x(m,i) is the element in position m,i of the measured matrix X. The operators ∠, Im⁡· and Re⁡· stand for the phase operator, the imaginary part of a complex number and the real part of a complex number, respectively.Since each channel is related to each line of the matrixX, in order to compute the phase shift between two channels, we compute the phase difference of two consecutive antennas(7)ωm,:=Φm,:-Φr,:∈C1×N,where Φ stands for the matrix containing the calculated phases by (6), r indicates the reference channel, and m varies from 1 to M. This reference channel can be randomly selected from 1 to M and is the input of the SDR that is used as a reference to compensate the phase imbalance from the other inputs. Since the vector ω(m,:) is the m-th row of matrix Ω, in case m=i, the i-th row of Ω is filled with zeros. Finally, since the phase difference may slightly vary for different samples in the same row of Ω due to the thermal noise, we compute the arithmetic mean of the elements of each row of Ω, obtaining the vector ω¯∈CM×1 and its m-th element is given by(8)ω¯m=1N∑n=1NΩm,n.Hence, in order to compensate the phase shift between two different channels, we compute the vectorc∈CM×1. The m-th element is given by(9)cm=e-jω¯m.Note that the compensation vectorc is computed only once for the system initialization. The calibrated outputs of the antenna array are given by the following expression:(10)Xc=diag⁡cX,where the operator diag· transforms its argument vector into the main diagonal of a diagonal matrix.


## 3.2. Hardware Assembling for the Drone Tracking Measurement Campaign

After the hardware has been calibrated, the next step is to assembly it in order to perform the measurements.The four-element omnidirectional antenna array is connected to the calibrated hardware composed of the FPGA motherboard and SDR daughterboard according to Figure2. Each antenna is dual band (from 2400 MHz to 2483.5 MHz and from 4900 MHz to 5875 MHz ) [36] and has linear polarization with 3.7 dBi of gain. The space Δ between two consecutive antennas is equal to 5.99 cm. The operational frequency f=2.5GHz is the maximum frequency that avoids aliasing.Figure 2
Proposed low cost antenna array based drone tracking device with 4 element ULA and a 4× 4 MIMO SDR.


## 3.3. Framework for DoA Estimation

Here, we first propose a sample selection approach for DoA estimation by automatic phase deviation detection. Then, we formulate a DoA estimation framework exploiting preprocessing techniques and model order selection schemes.Figure3 depicts the flowchart of the proposed signal processing solution for DoA estimation.Figure 3
Flowchart of the proposed solution for DoA estimation.As shown in Box 2 of the Figure3, the phase deviation correction proposed in Section 3.1 returns a matrix Xc that is used in the sample selection scheme in Section 3.3.1.


### 3.3.1. Sample Selection for DoA Estimation by Automatic Phase Deviation Detection

As exemplifiedXc in Figure 4, empirically we observed that the hardware causes phase deviations on the samples in random time instants. Therefore, we propose an approach to select the samples with phase deviations for the DoA estimation.Figure 4
Phase of the received data in each antenna. An example of phase deviation region is highlighted in the zoom area.Note that the phase compensation proposed in Section3.1 has been applied to the samples, whose matrix Φ containing the phases are depicted in Figure 4. Furthermore, note that there are significant deviations that can degrade the DoA estimation process. The main objective here is to remove these phase deviations.As shown in Figure5, such ripples can be better visualized by computing the phase difference in the time dimension according to the following expression:(11)γm,i=ϕm,i+1-ϕm,i2,where γ(m,i) is the value containing the quadratic difference of two consecutive time samples i and i+1 of the m-th channel. The ϕ(m,i), from (6), is the element in position (m,i) of the matrix Φ.Figure 5
Difference between consecutive samples from Figure4.Figure5 draws the Γ that contains the result of (11). By detecting the peaks, we can identify which samples should be removed. For this task, we apply the approach proposed in [37], which returns the green curve with the value of the threshold. Therefore, the samples whose phase differences are greater than the threshold are removed. The result after the samples removed is presented in the following equation:(12)Xcs=Xc:,1:N1∣Xc:,N2:N3∣⋯∣Xc:,NT-1:NT,where Xcs is the matrix with the selected samples. Note that NT=N and the values of Nt, for t=1,…,T, are found by comparing the phase difference values with the threshold in Figure 5.


### 3.3.2. DoA Estimation Framework

According to Figure3, the matrix Xcs given in (12) is used to improve the DoA estimation with preprocessing schemes. There are several DoA estimation schemes in the literature, such as beamforming approaches and subspace-based approaches. Examples of beamforming are Delay and Sum [38] and Capon [39], whereas examples of classical subspace-based approaches are MUSIC [23] and ESPRIT [40].The DoA estimation schemes assume that the model orderd is known. In practice, model order selection techniques should be applied to estimate the model order d, as depicted in Figure 3.In the flowchart of Figure3, we adopted the Exponential Fitting Test (EFT) [43, 44] as the model order selection scheme. The EFT has the deflation property that allows us to find suitable thresholds as a function of the Probability of False Alarm (Pfa). By exploiting the deflation property and by finding suitable thresholds, the EFT has been the only scheme in the literature to estimate d=1 in the presence of a strong LoS signal and d=0 in the only noise (no signal) measurements. We have compared several schemes in the literature such as Akaike Information Criterion (AIC) [45], Efficient Detection Criterion (EDC) [46], Minimum Description Length (MDL) [47], Stein’s unbiased risk estimate (SURE) [48], RADOI [49], ESTimation ERror (ESTER) [50], and Subspace-based Automatic Model Order Selection (SAMOS) [51]. The M-EFT [43, 44, 52] has also been suitable, but an even smaller Pfa was required to find the thresholds. The computation of the thresholds of the EFT requires an extremely low Pfa. The complexity of such a computation is prohibitive. Therefore, we propose in Appendix A an extrapolation algorithm to compute such thresholds. Note that our proposed extrapolation algorithm has been applied in [53–56], although no details are provided. The reason for extremely low Pfa may be related to the noise behaviour as shown in Appendix B. Note that the Ilmenau Package for Model Order Selection (IPM) [57] with MATLAB and Java implementation of the model order selection schemes can be found at the LASP homepage (https://lasp.unb.br/index.php/publications/softwares/).In order to further improve the accuracy of DoA estimation schemes, preprocessing schemes can be applied beforehand. We consider in this work the Vandermonde Invariation Technique (VIT) [58], Spatial smoothing (SPS) [59], and Forward Backward Averaging (FBA) [60, 61] as a preprocessing schemes. As depicted in Figure 3, after the preprocessing step, a matrixZ is returned and used by the DoA methods summarized in Table 1.Table 1
Selected state-of-the-art DoA estimation schemes.
Delay And Sum [38] P D S θ = w θ H R x x w θ w θ H w θ . ( 13 )   Capon [39] P C A P ( θ ) = 1 w θ H R x x - 1 w θ . ( 14 )    MUSIC [23] P M U S I C ( θ ) = 1 w H θ V n V n H w θ, ( 15 )   ESPRIT [40] Ψ = J 1 U s + J 2 U s ( 16 ) Ψ = E V D E Φ E H ,   w i t h   Φ = d i a g  [ ϕ 1 , … , ϕ d ] ( 17 ) μ i = ∠ ( ϕ i ) ,   i = 1 , … , d ( 18 )In Table1, the vectorw(θ) in (13), (14), and (15) vary according to the candidate values of θ. The value of θ that maximizes the expression in (13), (14), and (15) is the θ^1, since, in Section 2, the data model assumes that the LoS component faces no obstacles. Therefore, the component corresponding to the greatest power should be the same component with DoA equal to θ^1. Us∈CM×d is signal subspace, which is equal to the d eigenvectors corresponding to the d greatest eigenvalues, whereas Vn∈CM×M-d is a basis for the noise subspace, composed by the M-d eigenvectors associated with to the M-d smallest eigenvalues. In (17) and (18), Φ is the diagonal matrix that has the eigenvalues of Ψ. We compute all the spatial frequencies and the one whose component has the greatest power is the μ^1 that can be mapped to θ^1.


## 3.3.1. Sample Selection for DoA Estimation by Automatic Phase Deviation Detection

As exemplifiedXc in Figure 4, empirically we observed that the hardware causes phase deviations on the samples in random time instants. Therefore, we propose an approach to select the samples with phase deviations for the DoA estimation.Figure 4
Phase of the received data in each antenna. An example of phase deviation region is highlighted in the zoom area.Note that the phase compensation proposed in Section3.1 has been applied to the samples, whose matrix Φ containing the phases are depicted in Figure 4. Furthermore, note that there are significant deviations that can degrade the DoA estimation process. The main objective here is to remove these phase deviations.As shown in Figure5, such ripples can be better visualized by computing the phase difference in the time dimension according to the following expression:(11)γm,i=ϕm,i+1-ϕm,i2,where γ(m,i) is the value containing the quadratic difference of two consecutive time samples i and i+1 of the m-th channel. The ϕ(m,i), from (6), is the element in position (m,i) of the matrix Φ.Figure 5
Difference between consecutive samples from Figure4.Figure5 draws the Γ that contains the result of (11). By detecting the peaks, we can identify which samples should be removed. For this task, we apply the approach proposed in [37], which returns the green curve with the value of the threshold. Therefore, the samples whose phase differences are greater than the threshold are removed. The result after the samples removed is presented in the following equation:(12)Xcs=Xc:,1:N1∣Xc:,N2:N3∣⋯∣Xc:,NT-1:NT,where Xcs is the matrix with the selected samples. Note that NT=N and the values of Nt, for t=1,…,T, are found by comparing the phase difference values with the threshold in Figure 5.


## 3.3.2. DoA Estimation Framework

According to Figure3, the matrix Xcs given in (12) is used to improve the DoA estimation with preprocessing schemes. There are several DoA estimation schemes in the literature, such as beamforming approaches and subspace-based approaches. Examples of beamforming are Delay and Sum [38] and Capon [39], whereas examples of classical subspace-based approaches are MUSIC [23] and ESPRIT [40].The DoA estimation schemes assume that the model orderd is known. In practice, model order selection techniques should be applied to estimate the model order d, as depicted in Figure 3.In the flowchart of Figure3, we adopted the Exponential Fitting Test (EFT) [43, 44] as the model order selection scheme. The EFT has the deflation property that allows us to find suitable thresholds as a function of the Probability of False Alarm (Pfa). By exploiting the deflation property and by finding suitable thresholds, the EFT has been the only scheme in the literature to estimate d=1 in the presence of a strong LoS signal and d=0 in the only noise (no signal) measurements. We have compared several schemes in the literature such as Akaike Information Criterion (AIC) [45], Efficient Detection Criterion (EDC) [46], Minimum Description Length (MDL) [47], Stein’s unbiased risk estimate (SURE) [48], RADOI [49], ESTimation ERror (ESTER) [50], and Subspace-based Automatic Model Order Selection (SAMOS) [51]. The M-EFT [43, 44, 52] has also been suitable, but an even smaller Pfa was required to find the thresholds. The computation of the thresholds of the EFT requires an extremely low Pfa. The complexity of such a computation is prohibitive. Therefore, we propose in Appendix A an extrapolation algorithm to compute such thresholds. Note that our proposed extrapolation algorithm has been applied in [53–56], although no details are provided. The reason for extremely low Pfa may be related to the noise behaviour as shown in Appendix B. Note that the Ilmenau Package for Model Order Selection (IPM) [57] with MATLAB and Java implementation of the model order selection schemes can be found at the LASP homepage (https://lasp.unb.br/index.php/publications/softwares/).In order to further improve the accuracy of DoA estimation schemes, preprocessing schemes can be applied beforehand. We consider in this work the Vandermonde Invariation Technique (VIT) [58], Spatial smoothing (SPS) [59], and Forward Backward Averaging (FBA) [60, 61] as a preprocessing schemes. As depicted in Figure 3, after the preprocessing step, a matrixZ is returned and used by the DoA methods summarized in Table 1.Table 1
Selected state-of-the-art DoA estimation schemes.
Delay And Sum [38] P D S θ = w θ H R x x w θ w θ H w θ . ( 13 )   Capon [39] P C A P ( θ ) = 1 w θ H R x x - 1 w θ . ( 14 )    MUSIC [23] P M U S I C ( θ ) = 1 w H θ V n V n H w θ, ( 15 )   ESPRIT [40] Ψ = J 1 U s + J 2 U s ( 16 ) Ψ = E V D E Φ E H ,   w i t h   Φ = d i a g  [ ϕ 1 , … , ϕ d ] ( 17 ) μ i = ∠ ( ϕ i ) ,   i = 1 , … , d ( 18 )In Table1, the vectorw(θ) in (13), (14), and (15) vary according to the candidate values of θ. The value of θ that maximizes the expression in (13), (14), and (15) is the θ^1, since, in Section 2, the data model assumes that the LoS component faces no obstacles. Therefore, the component corresponding to the greatest power should be the same component with DoA equal to θ^1. Us∈CM×d is signal subspace, which is equal to the d eigenvectors corresponding to the d greatest eigenvalues, whereas Vn∈CM×M-d is a basis for the noise subspace, composed by the M-d eigenvectors associated with to the M-d smallest eigenvalues. In (17) and (18), Φ is the diagonal matrix that has the eigenvalues of Ψ. We compute all the spatial frequencies and the one whose component has the greatest power is the μ^1 that can be mapped to θ^1.


## 3.4. Tensor-Based DoA Estimation

In this subsection, a tensor factorization, namely, the PARAFAC decomposition, is applied. The PARAFAC decomposition generates three factor matrices from a received tensorX, whose structure is detailed in this subsection. One factor matrix corresponds to the estimate of the steering matrix A containing DoA information. In this subsection A is extracted from tensor X to estimate the DoA of the impinging signal.We first consider an unchanging sequence ofN symbols transmitted periodically. Such symbols can be found in a header or footer or even in the payload of a message. Alternatively, repeating sequences of symbols can be extracted from time periods when no data is being transmitted, but the carrier of the transmitter is active. An example is given as follows. Since the oscillators at the transmitter and at the receiver are never exactly the same, a small frequency deviation or constant phase change is observed at the receiver. At an MSK receiver, if the deviation is positive, a sequence of ones can be extracted and if the deviation is negative a sequence of zeros is observed.The symbols corresponding to the sequence from thei-th source is represented by the vector ci∈CN×1. Accordingly, we can build a received signal matrix19Xip=γipaiμiciT,∈CM×Nwhere p is the period corresponding to the sequence transmission. Generalizing for d signals we find20Xp=ADγpCT,∈CM×N,where Dγ(p)=diag(γ1(p)γ2(p)…γd(p))T and C=[c1c2…cd].ForP transmitted sequences, the symbols can be concatenated along the third dimension to form the received signal tensor21X=X1⊔3X2⊔3…⊔3XP∈CM×N×P,where ⊔3 represents concatenation along the third dimension. Since the slices of X can be written as 20, X has a PARAFAC structure and can be decomposed into three factor matrices A, C and Γ, where Γ contains the diagonals of Dγ(p) along its rows. To factorize X, we first rewrite it in three different matrix representations or unfoldings22X1=AC⋄ΓT,23X2=CΓ⋄AT,24X3=ΓA⋄CT,where ⋄ is the Khatri-Rao (column-wise Kroenecker) product.It is known that minimizing22, 23, and 24 in the least squares sense [62] leads to the following solutions:25A^=X1C⋄ΓT+,26Γ^=X2Γ⋄AT+,27C^=X3A⋄CT+.We consider the well-known Alternating Least Squares (ALS) algorithm to solve25, 26, and 27 in an iterative way. Since C is known, the estimation step 27 is skipped, and the ALS algorithm alternates between the estimations of A and Γ in a two-step approach [63].Once convergence is achieved, we useA^ to extract the DoA of the i-th source as follows: 28θ^i=argmaxθ⁡aiHθa^i


## 4. Experiments

In this section, we validate our proposed drone tracking device with measurement campaigns in an outdoor scenario. In Section4.1, the setup for the measurement campaign is described, whereas, in Section 4.2, we present the obtained results.


### 4.1. Experimental Setup

In Figure6, we depict the outdoor scenario used for the measurement campaigns. On the right-hand, we placed our drone tracking device, proposed in Section 3, as the receiver, whereas, on the left-hand, the transmitter is placed. The transmitter is a 2x2 MIMO SDR platform ASPR4 [64], with frequencies ranging from 50 MHz to 6.0 GHz, a channel bandwidth varying from 200 kHz to 56 MHz and a maximum power of 10 dBm at the output port.Figure 6
Photo taken during the measurement campaigns of the outdoor scenario pointing out the positions of the transmitter and of the antenna array based receiver.As shown in Figure7, the distance between the transmitter and the receiver is 48 m. Both transmitter and receiver are placed on tripods 115 cm above the ground. Note that the red “X” in Figure 7 is the location from where the photo in Figure 6 has been taken.Figure 7
Top view of the outdoor scenario for the measurement campaigns including the positions of the transmitter and of the antenna array based receiver.Both the transmitter and the receiver were set up using a MSK message signal, at2.48 GHz carrier frequency and 250 kbps of data rate. Before starting the experiment, the 2.48 GHz frequency was scanned and no noise source was detected. In order to verify that the receiver properly works and measures the Bit Error Rate (BER), we have to decode the signal. To this purpose, the transmitted package must be known and consists of pseudo random sequences of length 1024 bits and a header and footer with 16 bits each defined as 0xFFFF and 0x0000, respectively. Therefore, the total size of the package is 1056 bits. The transmitter uses both sampling frequency and bandwidth of 2MHz. At the receiver side, a 4 MHz sampling frequency and 4 MHz bandwidth are used. Each captured frame at the receiver has 5120 samples.As shown in Figure7, the transmitter is fixed and the receiver rotates from +90° to -90° in steps of 10°. According to Table 2, which presents the number of captured frame per DoA, at each 10° step, about 31 frames of size 5120 samples are captured.Table 2
Number of frames captured at each10° step of the measurement campaign.
DoA - 90 ° - 80 ° -70° - 60 ° - 50 ° - 40 ° - 30 ° - 20 ° - 10 ° 0 ° 10 ° 20 ° 30 ° 40 ° 50 ° 60 ° 70 ° 80 ° 90 ° Frames 31 29 29 30 33 34 32 29 24 30 32 24 30 33 26 31 43 41 31


### 4.2. Experimental Results

This subsection shows the performance of the DoA estimation schemes shown in Table1. During the measurement campaign, the achieved Bit Error Rate (BER) was 10-4.Figure8 shows the DoA calculated by rotating the base array over the time. Our proposed device works for a DoA ranging from -60° to 60°. The DS, Capon, MUSIC, ESPRIT, and Tensor combined with preprocessing schemes are shown in Figure 9(b).Figure 8
Comparison between the DoA estimation schemes by varying the DoA from+60° to -60° with steps of 10°.Figure 9
Comparison of the DoA estimation error between the DoA estimation schemes from Table1 by varying the angle from -90° to 90° with steps of 10°.
(a)
Without preprocessing (b)
With preprocessingIn Tables3 and 4, we present the Root Mean Square Error RMSE for the schemes in Table 1 with and without preprocessing schemes using the measurements from Figure 6, respectively. The equation for RMSE is given by29RMSEθ1=1Q∑q=1Qθ1q-θ^1q2where q is one realization of a total of Q realizations for each 10° stage of the measurement campaign. In other words, the acquired data at each 10° degree step is reshaped into small matrices. The computation of the DoA is individually performed for each of these matrices. In our experiment, we chose empirically matrices of size 4 × 1000. Since, as shown in Table 2, at each 10° step 31 frames of size 5120 samples are captured, in average there are approximately Q=158 matrices depending on the data reduction performed in Section 3.3.1. The variables θ and θ^ stand for the actual and the estimated DoA, respectively.Table 3
RMSE for the schemes in Table1 without the preprocessing using the measurements from Figure 8.
Algorithm DS CAP. MUS. ESP. TEN. RMSE 2.3 ° 1.7 ° 4.0 ° 2.8 ° 1.7°  Variance of the RMSE 4,3° 1.1 ° 11.5 ° 4.1 ° 4.4°Table 4
RMSE for the schemes in Table1 with the preprocessing using the measurements from Figure 8.
Algorithm DS CAP. MUS. ESP. TEN. RMSE 1,6° 1,7° 1,7° 1,9° 2,7°  Variance of the RMSE 2,3° 2,5° 2,4° 1,8° 23,7°Comparing Tables3 and 4 we can note that, except for the Tensor, the algorithms presented improvement in terms of RMSE after incorporating the preprocessing. CAPON and the Tensor increased the variance with the preprocessing. The smallest RMSE was achieved by the DS approach after preprocessing. Note that the ESPRIT assumes the shift invariance property, while both MUSIC and ESPRIT exploit the property of the orthogonality between signal and noise subspaces. Note that both assumptions are approximations and, therefore, their exploitation may cause additional errors.In Table5, we compare the DoA estimation results obtained with our framework and with the state-of-the-art approaches. Furthermore, we show that, even without involving anechoic chamber in our tests, precise results are obtained in comparison to references that implemented its tests in a nonreflexive environment. The ‘x’ means that the commercial solutions did not provide such information. Finally, we provide the information about which of the solutions has low cost.Table 5
Comparison between the obtained results and the DoA estimation in practical measurements results already available in the literature.
Drone detect solutions General DoA solutions Solutions Proposed [25] [27] [20] [21] [22] [24] [15]  DoA precision 1.9 ° 1 ° 1 - 3 ° 5 ° 1 ° 4 ° 0.5 - 2 ° 2 °  Uses anechoic chamber No x x No Yes Yes Yes Yes  Low cost Yes No No Yes Yes Yes Yes YesAs shown in Table6, the total cost of the proposed drone tracking solution is US$ 2.222, whereas the solutions in [25] and [27] cost US$ 226.000 and US$ 120.000, respectively. Therefore, our proposed off-the-shelf solution costs less than 2 % of the commercial solutions in [25, 27].Table 6
Price table for the hardware of the proposed drone tracker device.
Hardware Cost (US$) AD-FMCOMMS5-EBZ-ND [33] 1080 4×Dual-Band Antenna 2.4 & 5GHz [36] 11 ZYNQ 7000 Zc702 [35] 999 Power Divider 2.4-6GHz 30Watts RoHS IP67 [41] 109 5×cable 305mm HPP100 SMA [42] 23 Total 2222


## 4.1. Experimental Setup

In Figure6, we depict the outdoor scenario used for the measurement campaigns. On the right-hand, we placed our drone tracking device, proposed in Section 3, as the receiver, whereas, on the left-hand, the transmitter is placed. The transmitter is a 2x2 MIMO SDR platform ASPR4 [64], with frequencies ranging from 50 MHz to 6.0 GHz, a channel bandwidth varying from 200 kHz to 56 MHz and a maximum power of 10 dBm at the output port.Figure 6
Photo taken during the measurement campaigns of the outdoor scenario pointing out the positions of the transmitter and of the antenna array based receiver.As shown in Figure7, the distance between the transmitter and the receiver is 48 m. Both transmitter and receiver are placed on tripods 115 cm above the ground. Note that the red “X” in Figure 7 is the location from where the photo in Figure 6 has been taken.Figure 7
Top view of the outdoor scenario for the measurement campaigns including the positions of the transmitter and of the antenna array based receiver.Both the transmitter and the receiver were set up using a MSK message signal, at2.48 GHz carrier frequency and 250 kbps of data rate. Before starting the experiment, the 2.48 GHz frequency was scanned and no noise source was detected. In order to verify that the receiver properly works and measures the Bit Error Rate (BER), we have to decode the signal. To this purpose, the transmitted package must be known and consists of pseudo random sequences of length 1024 bits and a header and footer with 16 bits each defined as 0xFFFF and 0x0000, respectively. Therefore, the total size of the package is 1056 bits. The transmitter uses both sampling frequency and bandwidth of 2MHz. At the receiver side, a 4 MHz sampling frequency and 4 MHz bandwidth are used. Each captured frame at the receiver has 5120 samples.As shown in Figure7, the transmitter is fixed and the receiver rotates from +90° to -90° in steps of 10°. According to Table 2, which presents the number of captured frame per DoA, at each 10° step, about 31 frames of size 5120 samples are captured.Table 2
Number of frames captured at each10° step of the measurement campaign.
DoA - 90 ° - 80 ° -70° - 60 ° - 50 ° - 40 ° - 30 ° - 20 ° - 10 ° 0 ° 10 ° 20 ° 30 ° 40 ° 50 ° 60 ° 70 ° 80 ° 90 ° Frames 31 29 29 30 33 34 32 29 24 30 32 24 30 33 26 31 43 41 31


## 4.2. Experimental Results

This subsection shows the performance of the DoA estimation schemes shown in Table1. During the measurement campaign, the achieved Bit Error Rate (BER) was 10-4.Figure8 shows the DoA calculated by rotating the base array over the time. Our proposed device works for a DoA ranging from -60° to 60°. The DS, Capon, MUSIC, ESPRIT, and Tensor combined with preprocessing schemes are shown in Figure 9(b).Figure 8
Comparison between the DoA estimation schemes by varying the DoA from+60° to -60° with steps of 10°.Figure 9
Comparison of the DoA estimation error between the DoA estimation schemes from Table1 by varying the angle from -90° to 90° with steps of 10°.
(a)
Without preprocessing (b)
With preprocessingIn Tables3 and 4, we present the Root Mean Square Error RMSE for the schemes in Table 1 with and without preprocessing schemes using the measurements from Figure 6, respectively. The equation for RMSE is given by29RMSEθ1=1Q∑q=1Qθ1q-θ^1q2where q is one realization of a total of Q realizations for each 10° stage of the measurement campaign. In other words, the acquired data at each 10° degree step is reshaped into small matrices. The computation of the DoA is individually performed for each of these matrices. In our experiment, we chose empirically matrices of size 4 × 1000. Since, as shown in Table 2, at each 10° step 31 frames of size 5120 samples are captured, in average there are approximately Q=158 matrices depending on the data reduction performed in Section 3.3.1. The variables θ and θ^ stand for the actual and the estimated DoA, respectively.Table 3
RMSE for the schemes in Table1 without the preprocessing using the measurements from Figure 8.
Algorithm DS CAP. MUS. ESP. TEN. RMSE 2.3 ° 1.7 ° 4.0 ° 2.8 ° 1.7°  Variance of the RMSE 4,3° 1.1 ° 11.5 ° 4.1 ° 4.4°Table 4
RMSE for the schemes in Table1 with the preprocessing using the measurements from Figure 8.
Algorithm DS CAP. MUS. ESP. TEN. RMSE 1,6° 1,7° 1,7° 1,9° 2,7°  Variance of the RMSE 2,3° 2,5° 2,4° 1,8° 23,7°Comparing Tables3 and 4 we can note that, except for the Tensor, the algorithms presented improvement in terms of RMSE after incorporating the preprocessing. CAPON and the Tensor increased the variance with the preprocessing. The smallest RMSE was achieved by the DS approach after preprocessing. Note that the ESPRIT assumes the shift invariance property, while both MUSIC and ESPRIT exploit the property of the orthogonality between signal and noise subspaces. Note that both assumptions are approximations and, therefore, their exploitation may cause additional errors.In Table5, we compare the DoA estimation results obtained with our framework and with the state-of-the-art approaches. Furthermore, we show that, even without involving anechoic chamber in our tests, precise results are obtained in comparison to references that implemented its tests in a nonreflexive environment. The ‘x’ means that the commercial solutions did not provide such information. Finally, we provide the information about which of the solutions has low cost.Table 5
Comparison between the obtained results and the DoA estimation in practical measurements results already available in the literature.
Drone detect solutions General DoA solutions Solutions Proposed [25] [27] [20] [21] [22] [24] [15]  DoA precision 1.9 ° 1 ° 1 - 3 ° 5 ° 1 ° 4 ° 0.5 - 2 ° 2 °  Uses anechoic chamber No x x No Yes Yes Yes Yes  Low cost Yes No No Yes Yes Yes Yes YesAs shown in Table6, the total cost of the proposed drone tracking solution is US$ 2.222, whereas the solutions in [25] and [27] cost US$ 226.000 and US$ 120.000, respectively. Therefore, our proposed off-the-shelf solution costs less than 2 % of the commercial solutions in [25, 27].Table 6
Price table for the hardware of the proposed drone tracker device.
Hardware Cost (US$) AD-FMCOMMS5-EBZ-ND [33] 1080 4×Dual-Band Antenna 2.4 & 5GHz [36] 11 ZYNQ 7000 Zc702 [35] 999 Power Divider 2.4-6GHz 30Watts RoHS IP67 [41] 109 5×cable 305mm HPP100 SMA [42] 23 Total 2222


## 5. Conclusions

In this paper, we have proposed a low cost antenna array based drone tracking device for outdoor environments. The proposed solution is divided into hardware and software parts. The hardware part of the proposed device is based on off-the-shelf components such as an omnidirectional antenna array, a 4-channel SDR platform with carrier frequency ranging from 70 MHz to 6 GHz, a FPGA motherboard, and a laptop. The software part includes algorithms for calibration, model order selection (MOS), and DoA estimation, including specific preprocessing steps to increase the DoA accuracy. We have evaluated the performance of our proposed low cost solution in outdoor scenarios. Our measurement campaigns have shown that when the array is in the front fire position, i.e., with a DoA ranging from-60° to 60°, the maximum and the average DoA errors are 6° and 1,9°, respectively. Our proposed off-the-shelf solution costs less than 2 % of commercial solutions in [25, 27]. In order to further improve our analysis of the proposed system and our results, experiments in an anechoic chamber can be performed. Moreover, the performance of the proposed framework can be improved by incorporating interpolation schemes. Perspectives also include the adoption of a more realistic noise model to simplify the computation of the thresholds of the Exponential Fitting Test (EFT).



---

*Source: 5437908-2019-06-27.xml*