An Ultra-fast DOA Estimator with Circular Array Interferometer Using Lookup Table Method

The time-consuming phase ambiguity resolution makes the uniform circular array (UCA) interferometer not suitable for real-time direction-of-arrival (DOA) estimation. This paper introduces the lookup table (LUT) method to solve this problem. The key of the method is that we look up the ambiguity numbers instead of the eventual DOA from the table, and then the DOA is obtained by relatively small amount of calculation. This makes it possible that we are able to shrink the table size while maintain the DOA estimation accuracy. The table addresses cover all possible measured phase differences (PDs), which enables the method to be free of spatial scanning. Moreover, without adding frequency index to the lookup table, the estimator can realize wideband application. As an example, a field-programmable gate array (FPGA) based DOA estimator with the estimation time of 180 ns is presented, accompanied by the measured results. This method possesses the advantages of ultra-high speed, high accuracy and low memory usage.


Introduction
DOA estimation is to determine the bearing of a radiation source and it has been widely applied in the military and the civilian fields including radar, sonar and communication.Different estimation methods have been proposed in the literatures [1], [2].They can be generally grouped into two categories.
The first is to make use of the amplitude information of signal received by the antenna [3].The signal entering the main beam of an antenna has the maximum power output.It has simple structure and algorithm but suffers from low accuracy.The second is to take advantage of the phase information of the signal.The interferometer and the well-known subspace type algorithms (i.e., MUSIC [4], ESPRIT [5]) both belong to this category.The subspace type algorithms are able to achieve super-resolution and detect multiple signals simultaneously.However, they are not suitable for the realtime application because of their high computational complexity and they require accurate array calibration [6,7].The interferometer has been widely applied in practice.DOA can be easily derived from the phase difference (PD) between antennas.It has high accuracy, simple structure and low cost.
For interferometer, when incident angle becomes larger, the phase difference may be out of the range [−φ, φ], then the so-called phase ambiguity will emerge.The relation between true phase difference and measured phase difference is φ t is the true phase difference and φ m is the measured phase difference which is in the range [−φ, φ].k is called the ambiguity number which is an integer.Regarding the ambiguity resolution, the combination of long baseline and short baseline is introduced on the linear array (LA) [8].The short baseline with its PD being unambiguous is responsible for resolving the ambiguity over the long baseline and the long baseline is used to keep the accuracy.Nevertheless, two LAs are needed for two dimensional (2-D) DOA estimation.They occupy too much space.Additionally, LA has its shortcoming for wideband application.The short baseline has to be short enough to make its PD unambiguous at high frequency.
To satisfy this, the antenna size has to be reduced sometimes, which may lead to deteriorated performance of antenna at low frequency.So, the circular array (CA) seems to be a good choice and there is the clustering based method for the ambiguity resolution on CA [9].Regrettably, the clustering based method is very time-consuming.
To avoid ambiguity resolution, the lookup table method is introduced in the correlative interferometer [10].It compares the measured PDs with the theoretical PDs whose mapping relationships with the true DOAs are stored in the table, and the bearing is obtained when the correlation coefficient is at a maximum.However, for 2-D DOA estimation, the spatial scanning is time-consuming.To improve the accuracy, the spatial grids have to be much denser and to realize wideband application, the table also has to be built along the frequency grid.Those will significantly increase the table size.
In this paper, we propose a LUT based method to estimate the DOA with a UCA interferometer.The timeconsuming ambiguity resolution is accelerated by the LUT method.Configuring the truncated PDs as the address, we look up the ambiguity numbers from the table.Then the ambiguities are resolved.With a small amount of calculation we will be able to achieve the DOA.Precisely because we look up the ambiguity numbers, the table size can be made relatively small while the DOA estimation accuracy is maintained and only the table at highest frequency within the band is needed for wideband application.Furthermore, the table addresses covering all possible PDs enables the method to be free of spatial scanning, which contributes to realizing ultra-high speed of DOA estimation.These will be explained in detail in Sec. 3. Before hardware implementation, some numerical simulations are conducted to determine the bit width of the table address and verify the performance.Finally, a FPGA based DOA estimator and its measured results are presented.Compared to the hardware implementation without using the LUT method reported in [11], our method speeds up the estimation by four orders of magnitude.Moreover, we also give some reasonable approaches to further improve the estimation speed.

Conventional Method for DOA Estimation with UCA
Consider a single signal impinging on a five-element UCA of radius R. A three dimensional Cartesian coordinate system is established on the array.As shown in Fig. 1, the angle γ = 18 • is used to locate the first sensor and ω = 72 • is the angle between sensors.The vector K is used to denote the signal whose DOA is (θ, ϕ).θ is the elevation angle which is in the range [0 • , 90 • ] and ϕ is the azimuth angle which is in the range [0 Assuming a noiseless scenario and setting the origin as the phase reference point, we can represent the output of the pth sensor at time t as where p = 1,..., 5 and i = √ −1. a p is the gain of the pth sensor, s(t) is the signal and λ is the wavelength of signal.In terms of the estimation accuracy, the long baseline is preferred.So we measure the PD between the pth and qth sensor, i.e., φ p,q (q = p + 2. if q = 6, let q = 1; if q = 7, let q = 2).Their true value can be shown as φ p,q =Arg(x p (t)/x q (t)) = − (4πR/λ) cos(θ) sin(ω) sin(ϕ − γ − ωp) where Arg(•) returns the phase angle of a complex variable.
When there is no phase ambiguity, we can directly use one pair of measured PDs to estimate the DOA.Taking (φ 1,3 , φ 2,4 ) as an example, first, we calculate the sum and the difference of them.
where ε = ϕ − γ − 3ω/2.Then we construct a complex variable f 1 which is At last, DOA can be obtained as follows However, in most cases phase ambiguities exist.Then a clustering based method is introduced to solve the ambiguities [9].
First, according to the desired frequency range, the desired field of view (FOV), array size and (3), we are able to obtain the maximum value of the true PD and then determine the range of ambiguity number by (1).Here we assume it to be in the range [−m, m].So the total number of possible ambiguity number is 2m + 1.
It is obvious that in each set there must be one and only one element that corresponds to the true DOA.In noiseless scenario, these five elements are the same.In noisy scenario they have the highest degree of clustering compared to the other elements selected from the five sets respectively.So our goal is to pick out the five most clustered elements, one element from one set.That can be done like this.Set one set as a reference (e.g., the first set).For each element of the set, find the four shortest Euclidean distance to the other four sets respectively and calculate the sum of the four distances.Then we will get (2m + 1) 2 summations.The element in the reference set having the minimum summation corresponds to the unambiguous PD pair and can be adopted to estimate the DOA through (7).At the same time we can output the true ambiguity numbers.
It can be seen that the clustering based method needs nested loops (loop within a loop) to resolve the ambiguities.The time complexity of ambiguity resolution is O (2m + 1) 4 (M − 1) where M denotes the number of sensors.It is very time-consuming and not suitable for realtime application.So we figure out a LUT method to accelerate it.

System Architecture
Figure 2 is the block diagram of LUT based DOA estimator with the five-element UCA interferometer.For ambiguity resolution we need five measured PDs which are φ 1,3 , φ 2,4 , φ 3,5 , φ 4,1 and φ 5,2 .However, one of them is redundant.For example, φ 5,2 can be derived by where mod(•) is the modulo operator.So we only have to measure four PDs (i.e., φ 1,3 , φ 2,4 , φ 3,5 and φ 4,1 ) and each of PDs is represented in L-bit.This can save the system cost and since the PDs are configured as the address of the table, the table size is shrunk.From the table, we look up the ambiguity numbers k 1,3 and k 2,4 instead of looking up the eventual DOA directly.This can bring several benefits.Firstly, Different from the DOA value which needs to be represented in large bit width, ambiguity numbers are small integers.So they can be represented in small bit width.Secondly, as the ambiguity numbers are the discrete variables, they are not sensitive to the accuracy of the measured PDs.That means if the measured PDs are represented in small bit width (e.g., N-bit and N < L), we are still able to achieve the true ambiguity numbers.So the N-bit PDs are configured as the address.When looking up the table, we truncate the L-bit measured PDs to N-bit by reserving the high N-bit.Thirdly, it is easy to find that there is no need to resolve the ambiguities at the exact frequency of signal.If the ambiguities are resolved at λ r , λ in (6) will be replaced by λ r and ( 6) is turned into The coefficient λ r /λ has no impact on the clustering based ambiguity resolution method described in section 2. Theoretically, the LUT method can achieve unlimited frequency bandwidth without adding the frequency component to the table address.All the above contribute to shrinking the table size.

k k
Fig. 2. Block diagram of LUT based method for DOA estimation.
After k 1,3 and k 2,4 being obtained, they are sent to the calculation unit which exerts (1), ( 4), ( 5), ( 6) and ( 7) to output the DOA.According to (1), ( 4), ( 5) and ( 6), we need to know λ and the measured φ 1,3 and φ 2,4 .λ can be achieved by the instantaneous frequency measurement (DIFM) [12] module in the radio frequency (RF) receiver and φ 1,3 , φ 2,4 obtained by the fast Fourier transform (FFT) or via a phase detector [13] is represented in full bit width (L-bit).As long as they are accurate enough, the DOA estimation accuracy will be maintained.The computational cost in the calculation unit is low as we can pre-calculate the constant components in (6).

Table Building
Table building is a procedure that we list all possible table addresses and for each calculate its corresponding table element.Since each PD configured as part of the address is an N-bit number, it has 2 N possible values which are in the set containing elements from −π to π with the interval of 2π/(2N − 1).Then the table address that concatenates the four truncated PDs together is a 4N-bit number and has 2 4N possible combinations.For each combination we will calculate the ambiguity numbers k 1,3 and k 2,4 with the method described in Sec. 2. The computational cost is high.However, the table is always pre-calculated and the calculation can be accelerated by the parallel technique.The table element is the combination of k 1,3 and k 2,4 .If it is represented by P bits, the table size will be 2 4N P/2 23 MB.Furthermore, there is still a question that which frequency the table is built at.The table built at the lower frequency within the band is a little different from the table built at the higher frequency.From (9) we can find that lower frequency has a possibility to make the magnitude of f 1r larger than 1.This f 1r will be removed during ambiguity resolution.If the f 1r is related to the true ambiguity numbers, then the table will be incorrectly built.Therefore the table should be built at the frequency greater than or equal to the highest frequency within desired band.Here we choose the highest frequency.

Numerical Simulations
Obviously, we want the bit width of the table address to be small enough to shrink the table size.Paradoxically, too small bit width may lead to unsuccessful ambiguity resolution.So, we conduct numerical simulations to determine N and verify the performance.
First, we define the DOA estimation error as where (•) denotes the inner product.K t is the vector associated to true DOA and K e is the vector associated to estimated DOA. Subscript t and e represent true and estimated, respectively.Suppose the table is built at R/λ = 2.  3 (a).DC means directly using the method in Section 2 without the LUT method.It can be seen that small bit width of table address decrease the successful probability of ambiguity resolution and N = 6 is a good choice since its performance is the nearly the same as the DC.So the bit width of the table address is 24.With the table element represented by 8-bit, the table size is 16MB.If the ambiguities are successfully resolved, the DOA estimation accuracy of the proposed method will be the same as the DC.So we leave out the accuracy comparison.
Next, we show the performance of the LUT method under wideband scenario.The simulation conditions are the same as previous except that the PDs are acquired by 128point FFT.R/λ is varied from 0.05 to 2.6 and N = 6.The root mean square error (RMSE) are depicted in Fig. 3 (b).We can find that as frequency goes low, the RMSE becomes large.This is because the effective array aperture is small at low frequency.So, in practical the bandwidth cannot be too wide.Then we test the range of FOV of the proposed method.It is determined by the elevation angle and the higher the frequency is, the smaller the FOV will be.So we vary θ at the highest frequency, i.e., R/λ = 2.6.Other condition is the same as above except that SNR is fixed at 10dB.Since the FOV here is approximately a circular cone, we depict its cross section in Fig. 4.
Different colors represent different RMSEs of the DOA estimation and their relationship is indicated by the color bar.It is easy to observe that when the elevation is in the range [0 • , 30 • ] the source can be successfully resolved.The FOV is limited due to the limited value of the maximum ambiguity number (i.e., m = 2) and can be enlarged by increasing the value of m.In practice, m should be determined in the first place by the desired FOV and frequency rang.

Advantages and Disadvantage
The advantages of the propose method include ultrahigh speed, high accuracy and low memory usage.
From Sec. 3.1 we see that the table size is shrunk by looking up the ambiguity numbers instead of the eventual DOA from the table.The memory usage will be lower than that of the correlative interferometer, especially for high accuracy wideband application.So the proposed method has the advantage of low memory usage.From the first simulation above, we find that when we set N = 6 the successful probability of ambiguity resolution of the proposed method is nearly the same as DC.And in the proposed method the DOA estimation accuracy will be the same as the DC, if the ambiguities are successfully resolved.Since the interferometer has the inherent advantage of high accuracy, we draw the conclusion that the proposed method also has high accuracy.
In contrast, the accuracy of the correlative interferometer depends on the density of spatial grids.For some DOA estimation methods, we need to perform spatial scanning.For instance, the correlative interferometer [10] needs to calculate the correlative coefficient at every spatial grid and the MUSIC algorithm [4] needs to calculate the value of a cost function for every direction vector.The spatial scanning can be very time-consuming, especially for the 2-D DOA estimation.However, in the proposed method, the table addresses cover all possible measured PDs due to the table building method described in Sec.3.2.So, without spatial scanning, we can look up the corresponding ambiguity numbers from the table directly.
Then the DOA can be derived by relatively small amount of calculation.So the proposed method is able to realize ultra-high speed.
However, in general, this method cannot deal with multiple signals simultaneously, which is the inherent disadvantage of the interferometer.But we can extract other informa-tion to separate the signals, e.g., we can tackle the multiple signals with different frequencies by spectral decomposition.

Practical DOA Estimator Design and Measured Results
In this section, a practical implementation of wideband DOA estimator is presented.We choose XC5VSX95T which is the Xilinx Virtex-5 FPGA as the processor and make it work at the clock of 200 MHz.The table is calculated and saved by MATLAB at the condition of N = 6, m = 2 and R/λ = 2.6.Then the table elements are written into a flash memory with their corresponding addresses.The flash memory we adopt is M29W128GH from Micron whose density is 16 MB.The hardware design is shown in Fig. 5.As the table address is in 24-bit, only the high 6-bit of each 14-bit PD are reserved by the truncator.From the flash memory we get the 8-bit data containing the ambiguity numbers k 1,3 and k 2,4 .Then the ambiguity numbers are sent to the calculation unit along with the 14-bit φ 1,3 and φ 2,4 .In the calculation unit, we make full use of the parallel computation and the pipeline technique to accelerate the calculation.It has a four-stage pipeline.At the first stage, the FPGA runs two independent (1) parallelly to output the true PDs φ t1,3 and φ t2,4 , respectively.At the second stage, (4) and ( 5) are also executed parallelly.At the third stage, the constant components (e.g., (8πR/λ) sin(ω) sin(ω/2)) in ( 6) can be pre-calculated.At last, the DOA which is represented in Q10 format in radian is achieved at the forth stage.
The hardware design is verified by Chipscope Pro Analyzer.For simplicity, the four measured PDs and the frequency associated with the DOA of (30 • , 60 • ) are passed the FPGA VIO Console.The waveform captured by in Fig. 6.PD valid means FPGA the PDs.After 14 clock cycles, the flash memory outputs the ambiguity numbers and gives a k valid pulse.The time between the PD valid pulse and k valid pulse depends on the access time of the flash memory.Afterwards, FPGA executes the operations in the calculation unit to output the DOA with a DOA valid pulse.The measured DOA in the figure is (536/2 10 , 1073/2 10 ) radian (i.e., (29.99 • , 60.04 • )) which is approximately equal to the true DOA.There are 36 clock cycles between the PD valid pulse and DOA valid pulse.As the clock frequency being 200MHz, the estimation costs 180 ns.There is a digital signal processor (DSP) based hardware implementation of the DC reported in [11].It takes 7ms to run a estimation.That is to say our LUT based method improves the speed by 38889 times, i.e., four orders of magnitude.
In addition, the estimation speed is able to be further accelerated by two approaches.The first is to increase the clock frequency with a more powerful FPGA (e.g., Virtex-6, Virtex-7).The second is to shorten the time table lookup costs by replacing the flash memory with the random access memory (RAM) and the table elements will be transferred to the RAM during stage of the system initialization.At last, we practically test the DOA estimation performance of the estimator in a microwave anechoic chamber.The antenna array is mounted on a servo and a radiation source whose wavelength satisfies R/λ = 2.2 is placed in the far-field region of the antenna.In our system, the four PDs are obtained by eight phase detectors AD8302 and then sampled by the four-channel 14-bit A/D convector AD7865.λ is measured by the DIFM.Since the servo can only swing in the vertical and horizontal plane, we define the the angle α and β to determine the source direction.α is the angle between the z-axis and the projection of soruce direction onto x-z plane.β is the angle between the z-axis and the projection of soruce direction onto y-z plane.we have tanα = cosφtanθ and tanβ = sinφtanθ.The servo scans the spatial domain and the DOA estimation results are uploaded to the upper computer.Then we depict the estimation error in Fig. 7.It can be see the source is successfully resolved at different location in the FOV.The failed regions at the four corners are out of the FOV.It should be note that the FOV here is larger than that in Fig. 4 because of the lower frequency.

Conclusion
This paper introduces the LUT method to accelerate the DOA estimation with a UCA interferometer.Owing to looking up the ambiguity numbers, the table is made relatively small without decreasing the accuracy.The table addresses covering all possible PDs enables the method to be free of spatial scanning.Monte Carlo simulations and FPGA based hardware implementation have confirmed its ultra high speed, high accuracy, and low memory usage.This method can also be extended to other array geometry.Since the lookup table address do not contain frequency component, the estimator can easily adapt to the wideband signal by the spectral decomposition.Therefore, the proposed method will be very attractive for the engineering applications.
6.The true DOA is (30 • , 60 • ).The wavelength of the signal is λ and R/λ = 2 and ambiguity number is in the range [−2, 2].The PDs are acquired by 64-point FFT and L = 14.Signalto-noise ratio (SNR) varies from −5 dB to 5 dB.For each condition, 500 Monte Carlo trials are run and the trail with DOA estimation error smaller than 3 • is assumed to be one successful ambiguity resolution.The simulation results of N = 4, 5, 6 and the direct calculation (DC) are shown in Fig.