Visibility Domain Direction-of-Arrival Optimization for Arbitrary 2-D Arrays

A visibility-modelling based method is introduced to determine the elevation and azimuth angles for a signal received at a quasi-randomly arranged (2-D) array. The method is optimal for the number of antennas when the array is configured such that no antenna baselines are repeated. Simulation shows that our method has comparable accuracy to subspace-based acrlong doa algorithms like acrfull music and measured data shows that it is both more accurate and faster than visibility based methods like acrlong ra all-sky imaging.


I. INTRODUCTION
Determining the direction-of-arrival (DoA) of radio signals incident on an antenna array is useful in many applications including the localization of radio-frequency interference (RFI), remote sensing, spectrum management, and next-generation communications systems. The standard techniques for DoA estimation, like beamforming, Multiple Signal Classification (MUSIC), and Estimation of Signal Parameters via Rotational Invariance Technique (ESPRIT), appeared simultaneously with the development of synthesis imaging in the Radio Astronomy (RA) community. As its name implies, synthesis imaging is intended to generate images of the radio sky using antenna arrays, but the formulation upon which it is built can be applied to DoA estimation as well.
The above mentioned DoA methods attempt to estimate DoA by determining the linear transformation between the incident signals and those received by the antennas as represented by a steering matrix. Unlike these methods, the visibility-modelling approach casts the cross-correlations between signals received by the array elements as samples of the two-dimensional Fourier transform of the incident radiation field. Estimating the direction of arrival of a signal is then equivalent to inverting this transformation, often done The associate editor coordinating the review of this manuscript and approving it for publication was Chengpeng Hao . by interpolating samples onto a monotonic grid and taking the 2-D Fast Fourier Transform (FFT), but achieved here through least-squares model fitting.
In this paper we show that compared to other methods, our proposed method is less computationally intensive, more noise resilient, handles near-horizon signals better, and is more accurate with fewer antennas.
The paper is structured as follows: in section I-A we review existing DoA methods as well as radio astronomy synthesis imaging techniques, in section II we present our method and in section III we validate the performance of our method on both simulated and measured data collected using the Long Wavelength Array (LWA). The incident signal is narrowband for both simulated and measured results.

A. LITERATURE
DoA estimation methods that work for truly arbitrary array geometries are either computationally expensive, like MUSIC [1], or suffer from low angular resolution, like conventional and Capon beamforming [2]. Faster variants of these algorithms exist, but rely on specific array geometries as in root-MUSIC [3], [4]. Other well-known methods like ESPRIT ( [5]) are based on pre-defined array geometry and are still often simplified to decrease computational load [6].
Much of the literature on these methods describes proposed algorithms as computationally efficient and available for arbitrary geometries and then exclusively use linear or FIGURE 1. Random array geometry for the Long Wavelength Array Sevilleta Station in New Mexico, a radio astronomy telescope. rectangular arrays, often with missing elements justifying the claim of arbitrary geometry [7]. We do not consider a linear or rectangular array with missing or slightly-shifted elements arbitrary. There are authors who call these arrays with missing or unequally spaced elements non-uniform [8], a choice in terminology we agree with. We will not review methods dependent on specific array geometries any further as we are only interested in truly arbitrary arrays.
For a method to work on an arbitrary array it should work for a random array. Random arrays are necessitated by a need to observe distant sources as in radio astronomy synthesis imaging (discussed below). Every baseline vector between two antennas that is repeated will yield a duplicate observation providing no further information about the source. Thus effort is spent minimizing the number of repeating baseline vectors during the array design process. There are several methods used for building random arrays such as a random position generator with constraints on maximum distance in any direction and minimum distance between elements [9]. Then each randomly generated array configuration is assessed for beam size/efficiency, radiation pattern, aperture efficiency, etc. This is often balanced with the cost of laying coaxial cable to each element. An example of a random array designed using this method is shown in fig. 1. Another random array design method for sparse (minimum element separation of λ > 1) arrays uses Poisson Disk sampling to achieve a similar effect as the random position generator while reducing the need for trial and error [10].
Synthesis imaging enables the creation of all-sky images using antenna arrays such as high-resolution radio telescopes [11]. This method relies on correlating the signals between each pair of antennas to obtain samples of the so-called visibility function.
Visibility is a function of array baseline coordinates (V (u, v)) whereas incident radiation intensity is a function of the direction cosines l, m as defined in eq. (1). The direction cosines are mappings from the celestial dome (described by elevation and azimuth) onto a flat plane. This relationship is illustrated in fig. 2. Just as traversing between the time and frequency domains is central to much of signal processing, moving between the incident radiation intensity and visibility domains is key to synthesis imaging [11], [12]. The theory of synthesis imaging and the Fourier transform relationship between these domains is reviewed briefly in appendix VI. The inversion of this transform to generate images or estimate properties of the incident radiation field is achieved in all-sky imaging through convolutional resampling of visibility samples onto a regular grid followed by a FFT.
Note that while a baseline spacing of 0.5λ may be optimal for traditional beamforming, for synthesis imaging many independent baselines are needed to fill the u, v space. This is also shown in appendix VI. The number of baselines required changes depending on the application: wideband radio telescopes like the LWA use hundreds of stands yielding hundreds of thousands of unique baselines to do space science with high angular resolution, yet much smaller 4-element arrays are using synthesis imaging to study ionospheric transients and structure [13].
In order to avoid the artifacts introduced by gridding, radio astronomers have also developed advanced techniques to identify cosmological entities like galaxies in visibility domain data [14]. These techniques are computationally slow and premised on extracting useful information about the entities like source position, shape, and flux parameters by fitting detailed models to the data. In contrast, we simplify the extracted information to only source location and thus do not require highly detailed visibility modelling.

II. METHODS
To determine signal DoA a parameterized model of the signal in the visibility domain is generated and its difference from VOLUME 10, 2022 the correlator output is minimized to determine the best-fit visibility model parameters. These parameters correspond to the source's position in the sky and can then be converted into the signal's DoA.

A. CONSTRUCTING A VISIBILITY MODEL
We model the signal received at the array as a point source with intensity I 0 and location (l 0 , m 0 ). The intensity and visibility domain models are then

B. FITTING TO SAMPLED VISIBILITIES
Fitting the model shown in eq. (3) to the visibilities output from the correlator (each correlator output is an integration, which is equivalent to a single snapshot in the language of MUSIC literature), is cast as least-squares minimization problem. The cost function has two parameters, l 0 and m 0 , and its minimum with respect to these parameters determines the best approximation of the source's location. For each of K measured points coming out of the correlator (in which case there are K baselines in the array) the residual r k is the difference between that measured value and the model. The least squares objective function is then the sum of the squares of the residuals at each measured visibility.
To minimize this non-linear function we elect to use the Levenberg-Marquardt algorithm for its ubiquity in the literature [11]. fig. 3 shows a contour plot of the objective function that the Levenberg-Marquardt algorithm must traverse to find the minimum at (l 0.45, m 0.45) for a signal at (θ = 0, φ = π 4 ). The function has steep gradients leading to the minimum but a near-flat domain otherwise as shown in fig. 3a. This requires an initial estimate (l i , m i ) accurate to 0.4 = (l i − l 0 ) 2 + (m i − m 0 ) 2 so that the Levenberg-Marquardt has a gradient to traverse. This is not appropriate for a general direction-finding method, and we propose a solution in section II-C.

C. CONSTRUCTING A DOMAIN-WIDE VISIBILITY MODEL
We bypass the need for an initial estimate by showing that the shape of the model we use during fitting does not necessarily need to be a close match to the shape of the data for the purpose of direction finding. We bypass the need for an initial estimate by showing that the shape of the model we are fitting does not matter in the context of direction finding. Although our received signal is best modelled by a point source (eq. (2) and (3), for the purpose of DoA identification we only need the peak of the model to align with the peak of the radiation intensity in the l, m domain. To this end we reconstruct our model so that the cost function is only flat at the minimum across the entire optimization domain.
We identify a 2-D Gaussian as a good choice for a domain-wide model with a defined peak that the Levenberg-Marquardt function can easily traverse.
The width of the Gaussian is defined by the width parameter α such that the Full Width at Half Maximum (FWHM) is It is clear from fig. 3c that an initial estimate of (0, 0) is always sufficient to fit a signal to a wide (α = 1) Gaussian model since the cost function domain is never flat. Whether a point-source model is used or a wide Gaussian model, this method only works for a single source since the objective function is traversed to the minima nearest the initial estimate. There are visibility-domain deconvolution methods like the widely used CLEAN algorithm which effectively handle multiple sources by iteratively locating and subtracting out sources [15]. Our method could be extended to handle multiple sources by doing a similar operation set to CLEAN, namely: 1) finding the center of a source, 2) fitting 2-D Gaussian parameters to the source, 3) subtracting the source out of the visibility domain data, 4) repeating for the next located source until there are no more sources (by either knowing the number of expected sources or defining a 2-D Gaussian fit threshold for α > 1 meaning it has fit to the residual noise in the visibilities). Via this extension, no apriori knowledge of the number of sources is needed. In this sense our proposed algorithm parallels the multi-signal capability of MUSIC. While the computed MUSIC spectrum has a peak for each received signal, a threshold must be defined above which a peak is considered a signal. This is trivial if the number of received signals is known but not otherwise. Similarly, the above algorithm can be run a pre-defined number of times if the number of signals is known but if it is not a Gaussian fit width threshold can be chosen above which the next fitted peak is considered noise.

III. RESULTS
To present our proposed method we first simulate data and validate the method accuracy against frequency, noise, number of antennas, integration time, and DoA. In all cases the model fitting is done using a 2-D Gaussian with α = 1 as described in section II-C. The methods chosen for comparison are MUSIC and all-sky imaging. While root-MUSIC and other MUSIC variants are computationally less expensive, we elect to use the original MUSIC algorithm for comparison as it has no array-geometry dependency and none of the variants provide greater accuracy than it.
We note that the literature surrounding MUSIC uses the term ''snapshot'' where we use the term integration (as does the radio astronomy field). While MUSIC calculates DoA based on a number of snapshots, we always calculate it based on a single integration. The controlled parameter in our method is the integration time.
Each simulation has different parameters but in all cases the signal is an unmodulated carrier wave as described by eq.
(2) or eq. (6) depending on whether the simulation takes place in the intensity domain or the visibility domain.

A. SIMULATION
For the following validations an unmodulated carrier wave with a planar wavefront was simulated and its phase at each antenna was computed from the antenna positions. These time series were run through the correlator to obtain the visibilities. To do so an array geometry is required and we use the antenna positions for the Long Wavelength Array Sevilleta Station (LWA-SV) available via the LWA Software Library (LSL) [16]. The benefit to using the LWA-SV array geometry is that we can directly compare our simulated results with measured data collected using the array (section III-B). The LWA-SV random array geometry is shown in fig. 1. Accuracy is measured as the angle between the input vector and the predicted/calculated DoA vector. This can either be the cos −1 of the dot product between the two vectors, or equivalently the magnitude difference between their l, m plane projections. Table 1 lists the simulations carried out and the values of the varied parameter for each.

1) ACCURACY VERSUS FREQUENCY/OFFSET
LWA-SV is built for operation from 4MHz to 88MHz and in a noiseless simulation our method has zero error caused by frequency across this bandwidth.
Our method relies on apriori knowledge of the received signal frequency and in a similarly noiseless simulation there is zero error for an offset between signal frequency and model fitting frequency of up to 100 kHz. Since the LWA-SV has an bandwidth of 100 kHz 1 this means our method does not see any error attributable to frequency offset using this array geometry.

2) ACCURACY VERSUS INTEGRATION TIME
LWA-SV returns a frame of samples from each antenna every 0.00512 s. This makes any integration time greater than this which is a multiple of 0.00512 s is easy to simulate using the LSL. When simulated in a noiseless environment there is zero error in accuracy using integration times ranging from 30 s to 0.00512 s.

3) ACCURACY VERSUS NUMBER OF ANTENNAS
A primary reason we use the LWA to validate our method is the high antenna count covering a large area (≈ 100 m × 110m). While the LWA consists of 520 dualpolarized antennas, we only make use of the 247 single polarization antennas oriented in the north-south direction. We only use a single polarization because our method depends on the spatial distribution of antennas and including both polarizations recorded at each point in space does not change the number of baselines we have to work with. Further we are interested in the DoA and not the polarization properties of the signal.
To simulate different array configurations we drop antennas from the correlator input and process those remaining. fig. 4 shows the accuracy curve as the overall array radius is held constant but the minimum baseline is increased. To do this ten different arrays were generated to satisfy each minimum baseline spacing requirement and the result averaged. The signal was set to the phase center of the array (θ = 90 • ). This is shown for a 10 MHz signal where MUSIC fails at ≈ 1.35λ (≈ 6 antennas) while all-sky imaging and our method continue to work until the minimum of 3 antennas (the constant error term for all-sky imaging is caused by the grid cell located around the array phase center). This result is consistent with the very good array performance obtained with only 4 elements in [13]. Not shown are similar results at 5 MHz and 15 MHz where all-sky imaging and our method continue to work down to the minimum of 3 antennas while MUSIC fails at 5 antennas. Note that the average number of antennas is not an integer. Multiple array configurations (with differing performance characteristics and number of antennas) are possible at each minimum spacing and we mitigate these effects by averaging 10 arrays at each minimum spacing. There is extensive literature discussing synthesis array optimization in the context of RA [11].

4) ACCURACY VERSUS SIGNAL ELEVATION
To compare the accuracy of this method as the signal DoA moves across the sky we step a simulated source from the array phase center (θ = 90 • ) to the horizon (θ = 0 • ) and assess the accuracy at each step. For both MUSIC and allsky imaging, as the signal approaches the horizon the sides of the celestial dome project onto a progressively smaller area of the l, m plane ( fig. 2), causing both an increasing distance between data points (grid cell centers) and a DoA bias towards the phase center as shown in fig. 5. This bias is not present with the proposed method. The saw-toothed profile of the all-sky imaging and MUSIC curves are due to the quantized l, m plane.
This results in the proposed method providing a substantial improvement over both MUSIC and all-sky imaging for lowelevation angle (θ near 0) signals. Signals of terrestrial origin refracted by the ionosphere will often arrive at low elevation angles. Thus an accurate DoA estimate in this range of angles is important for ionospheric studies.

5) ACCURACY VERSUS NOISE
In order to add noise to the simulation, white noise was added to the time series at each antenna as described by eq. (8)  where a is the time series peak amplitude for the complex signal s(t) = ae j2πft , snr is the desired Signal-to-Noise Ratio (SNR) as a linear ratio, SNR dB is the desired SNR in decibel and the resulting σ I and σ Q are the noise powers added to the real and imaginary components of the time series respectively.
This allows us to show accuracy versus SNR per sample at the correlator input as shown in fig. 6, a valuable metric for those operating in a communications system context.
With a 0.05 s correlator integration time we see that the accuracy against noise relation only becomes relevant at signal SNRs less than −20 dB. At signal powers below −20 dB all-sky imaging outperform our method which in turn outperforms MUSIC. However for signals with SNR greater than −23 dB our method is more accurate. Our method is fractionally worse than MUSIC for SNRs greater than −28 dB. This simulation was done at the phase center of the array (θ = 90) where it is apparent that all-sky imaging has a residual error, another artifact of the quantized l, m plane.
Below −36 dB (at 11 • error) our method stops returning DoA results and starts returning an out-of-bounds error. This means the model fit to a location where either of l or m is greater than 1. This begins to occur for all-sky imaging below −39 dB and does not occur in this SNR range for MUSIC. We believe this is a benefit of our method versus MUSIC and all-sky imaging as low accuracy results are clearly indicated.

B. MEASURED DATA
A 30 W radio transmitter was used to generate Continuous Wave (CW) signals with known DoA's received at the LWA-SV. It was operated from Santa Fe, New Mexico under the experimental radio license WK2XTU. The bearing across the earths surface between the transmitter and LWA-SV gives an azimuth based on geometry of 27.65 • . In practice, the azimuth may deviate from the geometric value as it reflects off the ionosphere. Due to these deviations we don't have a ground truth when considering measured data, however the simulations done in section III-A give us confidence that the method is accurate. As well, it is reasonable to assume the signal should have an azimuth near the geometric bearing.
MUSIC was omitted while processing hours of measured data because it proved prohibitively slow. Both our method and all-sky imaging make use of the very efficient LSL correlator, unlike our relatively naive MUSIC implementation. As well, MUSIC requires calculation of a pseudo-spectrum over the grid which is computationally intensive. fig. 7 shows the result of all-sky imaging and the proposed method for several hours of data recorded at LWA-SV. Results show close agreement with the geometric azimuth. The difference between mean calculated azimuth and geometric azimuth ranged from 0.20 • to 1.63 • for all-sky imaging and from 0.19 • to 1.10 • for our proposed method. The average standard deviation for all-sky imaging and our method closely agree in both azimuth and elevation as shown in Table 2. The deviations from geometric azimuth can be attributed to ionospheric perturbations and channel effects.
The gaps in the data appear where the received signal has a SNR low enough to cause an out-of-bound error (see section III-A5) or during periods where the experimental callsign WX2XTU is transmitted in Morse code for identification purposes.
The visible decrease in DoA uncertainty at higher signal frequencies ( fig. 7c versus fig. 7a) can likely be explained by the ionospheric channel. The highest frequency shown (7 MHz) is near the maximum usable frequency reflected by the ionosphere and thus there is likely only one signal path. Lower frequency signals generally have multiple paths The lower plot shows out-of-bound errors indicating that our method prefers returning a failure rather than a low accuracy result. manifesting here as DoA uncertainty (''noise'') in fig. 7. fig. 8 shows a histogram of azimuth values for our proposed method from the three observations shown in fig. 7 along with the geometric azimuth. This approximately Gaussian distribution with a mean at the geometric azimuth shows again that despite channel effects our DoA method functions as expected. The histograms are the same approximate shape independent of integration time but will have a varying sample count.
Processing a total of 6.9 h (24832 1 s integrations) of data yielded a consistently accurate result and showed that when run using similar computational resources our method was ≈ 18% faster than all-sky imaging.

IV. DISCUSSION
The results in fig. 7 show gradual changes in azimuth and elevation up to about 10 • over time periods of tens of minutes. These results may be interpreted as changes in the virtual height assuming a specular reflection from a mirror ionosphere.
Results not shown using shorter integration times of 0.1 s 0.5 s show no finer structure in the ionospheric virtual height measurements than is obtained with the 1 s integration VOLUME 10, 2022 The ≈ 18% speed-up from our method versus all-sky imaging is a rough estimate based on average computationtime-per-integration and over many hours of data ≈ 18% becomes substantial. In the pursuit of real-time processing (where the time to calculate DoA is less than integration time), there are also more opportunities to optimize our method as opposed to optimizing all-sky imaging by exploring new fitting models and optimization algorithms.
We wish to highlight that our proposed method is more of a framework than a strictly defined algorithm. Our chosen minimization algorithm (Levenberg-Marquardt) which traverses the objective function can be replaced with other traversal schemes. Similarly, our use of a wide 2-D Gaussian as a signal model when paired with the Levenberg-Marquardt removes the need for an initial estimate of signal DoA but can be replaced with any other signal model as befits the application. Changing the minimization algorithm and signal model while applying our method will yield different accuracy, resolution and computational load/time.
We recognize that our use of visibility domain model fitting is not novel as radio astronomers regularly do more detailed variants of this. Likewise we recognize that there exist sub-space based methods like MUSIC for identifying signal DoA on arbitrary 2-D arrays. We believe our contribution is another entry into the list of methods available for truly arbitrary array signal DoA.
The code for this work was written in Python and is available in a Git repository at https://github.com/mistic-lab/lwatools which contains methods for processing the measured data received at LWA-SV as well as the simulation framework used in this paper.

V. CONCLUSION
We have shown our proposed method is more accurate than MUSIC and all-sky imaging for near-horizon narrowband signals, while being comparable for variations in SNR, array radius, frequency and frequency offset. The simulated result for low elevation signals is notable as our method provides a substantially more accurate DoA measurement near the horizon than MUSIC and all-sky imaging. Terrestrial signals refracted by the ionosphere travel around the earth at low incidence angles and our method can be used to accurately determine their source direction. We then show results from running both our method and all-sky imaging on several hours of a narrowband signal collected at LWA-SV. Processing of the real data showed no loss of accuracy down to the minimum integration time of 0.1 s. This shows the validity of our method outside of a simulated environment as the mean of the resulting azimuth matches the across-earth bearing from transmitter to the receiver.

VI. DERIVATION OF COMPLEX VISIBILITIES
Consider two antennas which are separated by a baseline vector b. We can describe the baseline vector as a function of signal wavelength λ Assuming a planar array we neglect theẑ component. Let E 1 (l, m, t) and E 2 (l, m, t) denote the incident scalar electric field at the two antennas arriving from the direction indicated by the unit vectorŝ = lx + mŷ + √ 1 − l 2 − m 2ẑ at time t. l and m, thex andŷ elements ofŝ, are called direction cosines (eq. (1)), and can be written as functions of elevation (θ measured from the horizon, 0 • , up to zenith, 90 • ) and azimuth (φ measured from the local meridian, 0 • , clockwise while looking down). This is illustrated in fig. 2.
The intensity of the field incident on the array is I (l, m, t) = E 1 (l, m, t)E * 1 (l, m, t) where · denotes expected value with respect to time. This value is estimated using the integral of the product of the signals over a finite integration time τ .
I (l, m, t) = E 1 (l, m, t), E * 1 (l, m, t) τ The separation of the two antennas means that the incident wave fromŝ travels a distance of b · s between the two antennas, inducing a time delay in the signal received by the second. E 2 (l, m) is therefore E 2 (l, m, t) = E 1 (l, m, t)e j 2π λ b· s (12) Since the two antennas are sensitive to radiation from the entire sky, the signals seen at their outputs are, respectively, X 1 (t) = E 1 (l, m, t) dl dm X 2 (t) = E 2 (l, m, t) dl dm (13) The visibility associated with this pair of antennas is the correlation of their outputs expressed as a function of baseline coordinates V (u, v, t) = X 1 (t)X 2 (t) * τ = E 1 (l, m, t)E * 1 (l, m, t) τ e j2π(ul+vm) dl dm = Î (l, m, t)e j2π(ul+vm) dl dm (14) Exchanging the order of the expected value and integration requires that the source is spatially incoherent, which is trivially true for a point source. Every baseline b k given by a pair of antennas results in a single sample V (u k , v k ) of the continuous visibility function. eq. (14) takes the form of a 2-D Fourier transform and the incident radiation intensity can be recovered with an inverse Fourier transform.
In order to generate an all-sky image of I (u, v), the inverse transform must be evaluated on a regular grid of u, v points. In general this is computationally expensive, although resampling the visibility samples onto an appropriate grid permits the use of the FFT and a significant speedup. Visibility modeling, described above, avoids inverting this transform altogether.
It is also crucial to understand that the output of the correlator is integrated over some time τ in order to estimate the expected value of the product in equation 14. Whereas each antenna may be sampled at a high rate f s , the visibilities are delivered at a lower rate 1 τ . In that time the correlated values are repeatedly summed. In radio astronomy each output is called an integration whereas in MUSIC and similar algorithms it is known as a snapshot.