A Flexible Design Strategy for Three-Element Non-Uniform Linear Arrays

This paper illustrates a flexible design strategy for a three-element non-uniform linear array (NULA) aimed at estimating the direction of arrival (DoA) of a source of interest. Thanks to the spatial diversity resulting from non-uniform sensor spacings, satisfactory DoA estimation accuracies can be achieved by employing a very limited number of receiving elements. This makes NULA configurations particularly attractive for low-cost passive location applications. To estimate the DoA of the source of interest, we resort to the maximum likelihood estimator, and the proposed design strategy is obtained by constraining the maximum pairwise error probability to control the errors occurring due to outliers. In fact, it is well known that the accuracy of the maximum likelihood estimator is often degraded by outliers, especially when the signal-to-noise power ratio does not belong to the so-called asymptotic region. The imposed constraint allows for the defining of an admissible region in which the array should be selected. This region can be further modified to incorporate practical design constraints concerning the antenna element size and the positioning accuracy. The best admissible array is then compared to the one obtained with a conventional NULA design approach, where only antenna spacings multiple of λ/2 are considered, showing improved performance, which is also confirmed by the experimental results.


Introduction
Several sensors aim at estimating the direction of arrival (DoA) of narrowband signals exploit linear antenna arrays with a small number of receiving elements. This choice allows for the minimizing of the number of receiving channels and the reducing of the overall system complexity, enabling the design of low-cost, lightweight, and compact sensors which are suitable for the mass-market.
Among the many applications that could benefit from a small number of receiving channels are low-cost radars for vehicular anti-collision applications-which estimate the direction of the target echo passive sonar systems-as well as passive radar (PR) systems [1], usually characterized by low-cost requirements compared to their active counterparts. Specifically, PRs based on different RF waveforms of opportunity have been widely studied. For example, considering the drone detection application, encouraging results have been obtained by parasitically exploiting a variety of RF waveforms, such as GSM [2], UMTS [3], LTE [4], DVB-T [5][6][7][8][9], and Wi-Fi signals [10,11].
Furthermore, low-cost sensors are particularly attractive in civil surveillance applications, e.g., PR sensors for the surveillance of private facilities [12], passive location sensors for the indoor monitoring of human activity [13], and sensors for the localization of survivors during the search and rescue operations after natural disasters. For all these applications, cost containment, limited complexity, and system lightness are crucial features. Therefore, the number of employed antenna elements and receiving channels must be kept as low as possible.

1.
Achieving a longer global array length, thus improving DoA estimation accuracy; 2.
Attaining compliance with technological constraints (e.g., the antenna element size might set a minimum inter-element distance).
Recently, significant efforts have been invested in the optimization of this simple, but useful, configuration. Unfortunately, it is well known that when estimating the DoA in an angular sector [−θ 0 , θ 0 ], a ULA is affected by deterministic angular ambiguities if the uniform spacing between adjacent antenna elements exceeds λ/2[1 + sin θ 0 ]. In particular, in reference [20], the binomial weight configuration is used, and both the element distance and the array factor are selected for a specific case.
In contrast, a non-uniform linear array (NULA) is only affected by statistical ambiguities, whose impact largely depends on the available signal-to-noise power ratio (SNR). The use of NULA configurations has been explored in [21,22].
Specifically, under very high SNR conditions, the NULA estimation accuracy reaches the Cramer-Rao bound (CRB), which essentially depends on the global array length. In this SNR region, known as asymptotic region, the estimation variance is inversely proportional to the SNR. However, below a specific SNR value, the presence of statistical ambiguities results in large estimation errors, also known as outliers, occurring with a non-negligible probability. Several authors addressed this issue in literature [23][24][25][26][27], exploring the use of different lower bounds, such as the Ziv-Zakai bound [27]. In this range of SNR values, known as the threshold region, the outlier probability monotonically grows as the SNR decreases, while the DoA estimation variance rapidly deviates from the CRB.
Therefore, a design strategy to select the optimal NULA configuration is particularly attractive, since the operational application of the previously introduced low-cost sensors typically faces limited SNR conditions. The design strategy proposed in [28] provides useful guidelines to select a proper array configuration. However, only antenna spacings in multiples of λ/2 are considered, which in turn leads to the following two weaknesses:

1.
Depending on the SNR value and on the angular sector of interest, there is no guarantee that using inter-element distances with λ/2 quantization provides the optimal solution; 2.
Depending on the physical size of the employed antennas, operating with halfwavelength quantization may strongly limit the available design solutions, making the technological constraint more severe.
Moreover, it is unlikely that the quantized element spacing will allow for the investigation of the impact of the element positioning tolerance due to practical installation and/or manufacturing.
Array layout design approaches based on the symmetrical number system were proposed in [29,30], using an optimum symmetrical number system (OSNS) and a robust symmetrical number system (RSNS) of a pairwise relatively prime number, respectively. However, these design strategies do not provide a control of the performance in terms of statistical ambiguity resolution. More interestingly, an array spacing design procedure was presented in [31] for three-and four-elements arrays, including a constraint on a maximum admissible phase error in the design procedure, which allows for the exclusion of arrays that provide DoA ambiguity if the phase error exceeds a preassigned value. This approach can be considered as a first step towards our proposed design approach. In fact, in the following, we present an array spacing design procedure for the three-elements case that includes a direct control over the statistics of the ambiguity resolution. We fully characterize the probability of ambiguity resolution and provide an array spacing design procedure that guarantees a preassigned probability.
Specifically, based on the theoretical characterization of the DoA estimation error carried out by the author of [32], the proposed design strategy provides an array configuration which is not subject to the quantization constraint on the inter-element distance. The findings in [32] allow us to devise a procedure that constrains the searched NULA array configurations to provide an outlier probability below a maximum tolerable value. Moreover, the prediction of the maximum likelihood estimator (MLE) performance provided therein, in both the threshold and asymptotic regions, enables the identification of the solution with minimum mean square error (MSE), even for limited values of the available SNR. This also allows us to include additional constraints related to practical array manufacturing. It is interesting to note that design approaches for arrays with more than three elements were proposed in [33,34] as evolutions or combinations of the three-elements array procedures. Therefore, the contribution of this paper is also expected to provide useful elements for the design of longer arrays.
The remainder of the paper is organized as follows. In Section 2, the signal model, DoA MLE, and theoretical characterization of the MSE are introduced. In Section 3, the admissible pairwise probability (APP) region is defined as the parameters region containing the arrays that satisfy the constraint on the maximum outlier probability for a chosen SNR. Based on the APP region, the design strategy is derived as a constrained optimization problem and is presented in Section 4. In Section 5, a numerical analysis is carried out to assess the performance of various NULA configurations obtained with the proposed strategy. A comparison with the design strategy proposed in [27] is also carried out to show the limitations of the quantization constraint. Section 6 shows the benefits of the proposed strategy over the conventional approach when adding technological constraints, which further limits the available array configurations. Finally, in Section 7, the advantages shown for the simulated data are tested against an experimental dataset collected in the Wi-Fi band. Concluding remarks are reported in Section 8.

Array Geometry and Signal Model
Consider an N-element linear array receiving a narrowband signal from a single source with DoA θ s , measured relative to the array boresight. The displacements d n , n = 0, . . . , N − 1 of the element positions compared to the first element are collected into vector z = [0 d 1 . . . d N−1 ]. The described array geometry is sketched in Figure 1.
Let s(u) be the steering vector, defined as where λ is the wavelength, and u = sin θ is the steering direction. After I/Q demodulation to baseband, the N received signals are sampled and digitized. The snapshot of the array, collected at a specific time t 0 , is represented by the N × 1 complex vector x. Since the receiving elements are affected by thermal noise, the complex array output after downconversion, filtering, and sampling can be arranged into the N-dimensional column vector  Let ( ) be the steering vector, defined as where is the wavelength, and = sin is the steering direction. After I/Q demodulation to baseband, the received signals are sampled and digitized. The snapshot of the array, collected at a specific time , is represented by the 1 complex vector . Since the receiving elements are affected by thermal noise, the complex array output after down-conversion, filtering, and sampling can be arranged into the -dimensional column vector Assuming a deterministic signal model as in [32,35]: • is the generic sample of the complex baseband source signal and is modeled as a deterministic, but unknown, parameter. • ( ) is the target steering vector, obtained by Equation (1) along the direction of the source, namely = sin( ). • is an -dimensional column vector collecting the noise samples at the receiver.
Based on the pattern of the employed antenna elements and the possible presence of additional electromagnetic shielding, we assume that the antenna array can only receive signals with DoA in the sector = [− , ], being = sin( ). This angular sector is represented by the green region in Figure 1. Noise samples at the antenna elements are assumed to be distributed according to an independent, zero-mean complex Gaussian random variable, with the same mean square value , which is also independent of the source signal.

Maximum Likelihood DoA Estimation
To achieve the best estimation of the source DoA, we resort to the MLE. Based on the model in Equation (2), this is obtained by looking for the value of that maximizes the likelihood function Assuming a deterministic signal model as in [32,35]:

•
A is the generic sample of the complex baseband source signal and is modeled as a deterministic, but unknown, parameter. • S(u s ) is the target steering vector, obtained by Equation (1) along the direction of the source, namely u s = sin(θ s ). • n is an N-dimensional column vector collecting the noise samples at the receiver.
Based on the pattern of the employed antenna elements and the possible presence of additional electromagnetic shielding, we assume that the antenna array can only receive signals with DoA in the sector U max = [−u max , u max ], being u max = sin(θ max ). This angular sector is represented by the green region in Figure 1. Noise samples at the N antenna elements are assumed to be distributed according to an independent, zero-mean complex Gaussian random variable, with the same mean square value σ 2 n , which is also independent of the source signal.

Maximum Likelihood DoA Estimation
To achieve the best estimation of the source DoA, we resort to the MLE. Based on the model in Equation (2), this is obtained by looking for the value of u that maximizes the likelihood function Looking for the maximum in the region u ∈ U max , the ML estimate is given by: The estimateθ s of the DoA θ s is then obtained by inverting the relationship u s = sin(θ s ). Depending on the application at hand, the sources of interest may be only those with a DoA that lies inside the angular sector defined by U int = [−u int , u int ], being u int = sin(θ int ), possibly narrower than U max (namely θ int ≤ θ max ), as represented by the red region in Figure 1. Therefore, estimates that fall outside U int can be discarded as not being of interest. This allows for the removal of potential outliers of the ML estimation that correspond to largely displaced DoA. It is further noticed here that there might be a specific interest to only optimize the estimation performance inside an even narrower angular sector U opt represented by the yellow region in Figure 1, which will be discussed in the following sections of this paper.

Performance Evaluation and MSE Approximation
The design strategy presented in the following sections of this paper aims at selecting array configurations with the best possible performance. A basic building block for this procedure is provided by a robust performance prediction of the DoA estimation accuracy, evaluated in terms of MSE. A vast amount of technical literature has addressed the characterization of this MSE, which depends both on the array configuration and on the signal parameters. Figure 2 qualitatively shows the behavior of the estimation error: the continuous blue line represents the MSE, while the dashed blue line represents the CRB. gion in Figure 1. Therefore, estimates that fall outside can be discarded as not being of interest. This allows for the removal of potential outliers of the ML estimation that correspond to largely displaced DoA. It is further noticed here that there might be a specific interest to only optimize the estimation performance inside an even narrower angular sector represented by the yellow region in Figure 1, which will be discussed in the following sections of this paper.

Performance Evaluation and MSE Approximation
The design strategy presented in the following sections of this paper aims at selecting array configurations with the best possible performance. A basic building block for this procedure is provided by a robust performance prediction of the DoA estimation accuracy, evaluated in terms of MSE. A vast amount of technical literature has addressed the characterization of this MSE, which depends both on the array configuration and on the signal parameters. Figure 2 qualitatively shows the behavior of the estimation error: the continuous blue line represents the MSE, while the dashed blue line represents the . As shown in Figure 2, three different operative regions can be identified, which can be explained as follows. In the absence of noise, we have an ideal value of the likelihood function As shown in Figure 2, three different operative regions can be identified, which can be explained as follows. In the absence of noise, we have an ideal value of the likelihood function This provides the array beampattern for a source coming from u s , which always exhibits a maximum in the true target DoA, and has a mainlobe width that solely depends on the effective array length. For N > 2, L(u) is characterized by sidelobes, whose positions u m , m = 1, 2, . . ., levels, and widths depend on the array geometry (i.e., on the inter-element distances in z).
When the SNR is very high, there is a negligible probability that the noise values added to the signal source snapshot provide a peak of V(u) in the direction u m of a sidelobe of L(u), namely outside its main lobe, rather than in the direction of the true DoA, u s . In other words, in such conditions, thermal noise causes only slight variations in the peak location inside the main lobe, while the outlier probability is negligible. Therefore, the MSE of the ML DoA estimate accurately attains the CRB, derived as in [36][37][38][39] The CRB in Equation (6) only depends on the width of the main lobe, namely on the global array length, and is inversely proportional to the SNR. This region is referred to as the asymptotic region.
As the SNR decreases, the MSE rapidly deviates from its asymptotic behavior. Moreover, the probability that thermal noise results in a maximum of V(u) in correspondence of a sidelobe of L(u) becomes non-negligible. This leads to a DoA estimate well outside the main lobe of the likelihood function (namely an outlier), which results in a sensibly larger MSE. Therefore, the expression of the CRB can no longer be used as an approximation of the MSE. This threshold region is of major interest for practical applications, since in most cases, low-cost sensors need to operate with limited SNR values.
Finally, for very low SNR values, the probability that the effect of thermal noise provides a peak of V(u) in the sidelobes of L(u) becomes so high that the estimation process is undermined by the presence of several large errors. Under such conditions, the DoA estimates are uniformly distributed across the entire parameter space, and the MLE does not provide any reliable information on the true source DoA. This region is referred to as the no-information region.
Based on the discussion above, the authors of [32,35], derived an approximation of the MSE, denoted as ζ, achieved with the MLE in both the threshold and asymptotic regions: M represents the number of sidelobe peaks in the angular sector U int , where the estimates are considered of interest and are retained; • P m = P m (SNR, d 1 , . . . , d N−1 , u s ) represents the mth pairwise error probability, i.e., the probability that the mth sidelobe peak in position u m is chosen by the MLE search instead of the main lobe peak.
Notice that the union bound approximation derived in [40] and used in [32] ignores the events of two or more sidelobe peaks simultaneously being higher than the main lobe peak. Based on this, the probability of outlier P out can be approximated as Under a deterministic model assumption for the signal amplitude, a closed form expression for the pairwise error probability was also derived in [32] based on the discussion in [41], obtaining: where Q(a, b) is the Marcum Q-function, I 0 (·) is the modified Bessel function of the first kind and order 0, and its arguments depend on the array beampattern evaluated at the m-th sidelobe peak, located in u m = sin(θ m ): In the remainder of the paper, we build upon this theoretical performance characterization to derive a flexible design strategy for a three-element NULA.

Admissible Pairwise Probability Region
From the approximation in Equation (7), we notice that the values of P m act as scaling factors for the large error contributions (u m − u s ) 2 that add to the CRB to provide the global MSE. Therefore, the array configuration achieving the minimum MSE certainly lies in a region where the pairwise probabilities P m have negligible values. Consequently, an important step toward the array optimization consists of the identification of the region of the array parameters, namely d n , n = 1, . . . , N − 1, where all P m values are below an assigned threshold (i.e., the maximum acceptable value). This region will be referred to in the following as the APP region. To identify the APP region, we notice that the mth pairwise error probability in Equation (9) is a monotonic decreasing function of both the SNR and the array beampattern L(u m ) evaluated at the mth sidelobe peak location u m . Therefore, the array configuration achieving the minimum MSE depends on the specific SNR conditions. To carry out the optimization, we first select a fixed value for the SNR, i.e., SNR 0 . With this choice, the constraint on P m ≤ P max is equivalent to setting an upper bound on the value of L(u), namely L(u) ≤ L max , where L max = L max (SNR 0 , P max ).
By using the steering vector in Equation (1), the likelihood function L(u) in Equation (5) can be expressed in terms of the array parameters d n , n = 1, . . . , N − 1: where While this constraint is valid for any value of N, we are especially interested to the case of N = 3, which is appealing when addressing the design of low-cost, lightweight sensors, as described in the Section 1. For N = 3 array elements, the constraint in Equation (11) can be simplified as: which is a function of two parameters and can be easily represented in a plane. Figure 3 shows a contour plot representing the function L(∆ 1 , ∆ 2 ) on a ∆ 1 , ∆ 2 axis system.
which is a function of two parameters and can be easily represented in a plane. Figure 3 shows a contour plot representing the function (Δ , Δ ) on a Δ , Δ axis system. Notice that (Δ , Δ ) is a periodic function of both Δ and Δ with a unity period, i.e., (Δ − , Δ − ℎ) = (Δ , Δ ) for any ℎ, ∈ ℤ. For future reference, we define its single period as This is represented in the red box in Figure 3, where the numeric labels refer to the contour lines drawn for different possible values of (with ≤ = 9). Using the Notice that L(∆ 1 , ∆ 2 ) is a periodic function of both ∆ 1 and ∆ 2 with a unity period, i.e., L(∆ 1 − k, ∆ 2 − h) = L(∆ 1 , ∆ 2 ) for any h, k ∈ Z. For future reference, we define its single period L 0 as This is represented in the red box in Figure 3, where the numeric labels refer to the contour lines drawn for different possible values of L max (with L max ≤ N 2 = 9). Using the definition of L 0 (∆ 1 , ∆ 2 ), we can rewrite L(∆ 1 , ∆ 2 ) as so that the replica for k, h = 0 represents the main lobe of the ideal likelihood function, while the replicas for k, h = 0 represent its grating lobes.
Our goal now is to exploit the constraint on the likelihood function L(∆ 1 , ∆ 2 ) ≤ L max to identify the subset of array configurations satisfying the constraint on the maximum pairwise error probability, under the assumed condition of SNR = SNR 0 .
To achieve this, we first represent the beampattern L(u) of a specific array configuration z = [0, d 1 , d 2 ] in the plane (∆ 1 , ∆ 2 ). From Equation (12), we notice that ∆ 1 and ∆ 2 depend both on the array design parameters d 1 and d 2 and on the difference (u − u s ). However, the ratio between ∆ 2 and ∆ 1 is independent of (u − u s ), and it is equal to the ratio between d 2 and d 1 , . Therefore, as the steering direction u changes, the locus of points corresponding to the array beampattern L(u) of a specific 3-element array z = [0, d 1 , αd 1 ] can be represented in the plane (∆ 1 , ∆ 2 ) as a linear segment NP with slope α. For a specific array configuration, the endpoints N and P of the segment NP only depend on the maximum and minimum values assumed by the difference (u − u s ). By recalling that the array can only receive signals with DoA u s in the sector U max = [−u max , u max ], and that the we are only interested in estimating the DoA of the sources located in the sector U int = [−u int , u int ], we have u s ∈ U max and u ∈ U int . Hence, the difference u − u s is bounded as where µ max u max + u int . For a given µ max , the coordinates of the endpoints N and P for the array z in the plane (∆ 1 , ∆ 2 ) can be derived by inverting Equation (12) Given the biunivocal relationship between z and NP, determining whether an array configuration is admissible or not is now straightforward. Specifically, an array z can be considered admissible only if the corresponding segment in the (∆ 1 , ∆ 2 ) plane does not intersect the contour L 0 (∆ 1 − k, ∆ 2 − h) = L max for any k, h = 0. This condition guarantees that the sidelobes of the array beampattern L(u) are lower than L max , which in turn allows for satisfying the constraint on the maximum pairwise error probability. The following example allows for the visualization of the relationship between L(u) and NP, and further clarifies the notion of admissibility.
Let us assume an SNR 0 = 15 dB and a maximum acceptable pairwise probability P max = 10 −5 . By inverting Equation (9), we obtain L max ≈ 5.8612. The corresponding contour plot of the L(∆ 1 , ∆ 2 ) function is represented in Figure 4a. As is apparent, the region external to the contour lines, namely the green area, represents the set of (∆ 1 , ∆ 2 ) values that satisfy the constraint L(∆ 1 , ∆ 2 ) ≤ L max .    .
As explained previously, these array configurations can be represented in the plane (∆ 1 , ∆ 2 ) as segments lying on the line ∆ 2 = 2.1∆ 1 , and the extrema of these segments depend on the value of µ. Assuming µ max = 2 (i.e., u max = u int = 1), we evaluate the endpoints N and P through Equation (17), and we represent the four considered configurations in Figure 4a as colored segments NP to verify their admissibility.
As visible in Figure 4b, for ρ = 1, configurations z 1 and z 2 are admissible, since the segments representing their beampatterns do not intersect any replica of the likelihood function (except for the replica k = h = 0, which is the likelihood function main lobe). Conversely, configurations z 3 and z 4 are not admissible, since the segments representing their beampatterns do intersect the contour L 0 ( and for (k, h) = (−1, −2). The four array beampatterns are also represented in Figure 5 as a function of (u − u s ). Comparing these beampatterns, we notice that the array configurations denoted as z 1 and z 2 achieve the lowest sidelobes, while the others have high sidelobes that prevent fulfilling the requirement on L max .   Notice how the beampatterns of the array configurations characterized by the same α value are merely stretched versions of the same likelihood function. Specifically, the slope α = d 2 /d 1 defines the shape of the likelihood function, while the value µ determines the abscissa d 1 λ µ max of the ∆ 1 axis in which the likelihood function is cut to obtain the array beampattern.
As mentioned in Section 1, a straightforward approach to measure the robustness of a given array configuration to estimation ambiguities was proposed in [31], based on the maximum admissible phase error. Therein, the authors introduced the synthetic parameter min S pq , defined as the minimum distance between the folded replicas of the line φ 2 = αφ 1 , with the wrapped phase within the angular sector Φ = [0 • , 360 • ]. The minimum distance min(S pq ) of an array provides its maximum admissible phase error and is therefore a measure of the array's immunity to ambiguities. By assigning a minimum acceptable distance, the parameter min S pq can be used to determine the admissibility of an array.
The approach in [31] is inherently different from the one proposed in this work, since it capitalizes on the phase ambiguity of the involved baselines rather than on the statistical formulation of the MLE. However, the two approaches generally provide consistent results. To carry out a comparison, we used the algorithm described in [31] to evaluate min S pq for the array configurations in Figure 4a, with size scaled by ρ = 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, and reported their values, together with the corresponding maximum pairwise error probability max m P m , in Table 1 . Table 1, requiring that P m ≤ P max , ∀m is equivalent to requiring that min S pq > S min . However, while max m P m ranges from 10 −20 to 0.2, taking several intermediate values for different array configurations, min S pq only takes two values, namely 170.25 for more robust configurations, and 15.48 for less robust examples. This suggests that constraining min S pq provides an approximate admissibility criterion, while the proposed approach enables a more accurate assessment of array robustness to outliers. As a consequence, the two approaches provide consistent results for all array configurations, showing a well-defined behavior with respect to the presence of outliers. In contrast, they might yield conflicting indications for edge case configurations. Table 1. Minimum distance min S pq for the array configurations shown in Figure 4a and the corresponding maximum pairwise error probability max m P m . For example, from Table 1, we notice that the constraint P m ≤ P max = 10 −5 , ∀m is almost equivalent to the constraint min S pq > S min = 15.48. With these parameters, the two approaches agree on the admissibility of both configurations z 1 and z 2 for any ρ value. However, when ρ = 0.7, the criterion min S pq > 15.48 suggests that all the included array configurations should be considered admissible. Conversely, we notice that configuration z 4 does not satisfy the constraint P m ≤ 10 −5 , ∀m, being characterized by max m P m = 2 · 10 −5 , and thus it would be discarded using the proposed approach. Similar considerations apply to configuration z 3 when ρ = 1.

As visible in
Clearly, different values of SNR and P max can lead to different edge cases on these arrays. Furthermore, additional edge cases can be found, considering array configurations with different slopes α. However, the example above shows that the ML-derived proposed approach enables a more precise assessment of the array robustness to outliers, allowing us to better identify the boundaries of the APP region where the best constrained configurations in terms of MSE typically lie. This is shown in the following.
The green region in Figure 4b identifies the subset of abscissa values d 1 λ µ max , for which a cut of the likelihood function provides an admissible array beampattern, namely L(u), such that L(u) ≤ L max . The extrema of the green patch are denoted as ±∆ max 1 . This abscissa value defines the longest admissible array with α = 2.1, so that every array configuration is admissible. Notice that the value of ∆ max 1 depends on the chosen slope α. Hence, if we determine the value ∆ max 1 (α) of maximum abscissa for every possible slope α, we can identify the region in the plane (∆ 1 , ∆ 2 ) which contains every array configuration that satisfies L(u) ≤ L max , namely the APP region.
In principle, the value ∆ max 1 (α) can be derived as the intersection between the contours L 0 (∆ 1 − k, ∆ 2 − h) = L max and a line through the origin with slope α. Thus, the APP region can be derived in the plane (∆ 1 , ∆ 2 ) before setting a specific value for µ max .
Once the APP region has been identified, the value of µ max is required to relate the coordinates (∆ 1 , ∆ 2 ) to the inter-element distances (d 1 , d 2 ). In other words, as visible from Equation (17), the definition of the endpoints of the segment NP associated to an array z requires both the knowledge of (d 1 , d 2 ) and of µ max . Therefore, the admissibility of a specific array z depends not only on the constraint on the pairwise error probability, but also on how the sectors U max and U a are chosen.
For every slope α, the analytical value of the optimal ∆ max 1 (α) can be identified as follows. First, we notice that L 0 (∆ 1 , is symmetrical with respect to the two bisectors of the first and the second quadrants of the plane (∆ 1 , ∆ 2 ).
By defining the new set of coordinates γ and δ, such that: and applying the transformation, we obtain which can be easily solved for δ, obtaining: Using Equation (21) in Equation (19) allows us to express the points on the contour of L 0 (∆ 1 , ∆ 2 ) defined by L max , in terms of γ, as: Since the segment representing an array lies on a line through the origin, the slope of the line crossing the contour of the replica (k, h) is given by the ratio α k,h (γ) = ∆ 2 (γ)+h where α + k,h (γ) refers to the slope obtained from the positive solutions, and α − k,h (γ) refers to the slope obtained from the negative ones. As γ varies in the range (−2π, 2π), the results in Equation (23) provide the slopes of the lines crossing the replica (k, h) of the contour L 0 = L max and allow us to identify the maximum and minimum slope of its crossing line By exploiting the previous analytical re41sults, together with the observation that the lower order replicas block the limiting effect of the higher order replicas, we propose the following procedure to identify the APP region in the (∆ 1 , ∆ 2 ) plane. First, we identify the ∆ max 1 (α) relative to the slopes crossing the lowest order replica, then we progressively evaluate the ∆ max 1 (α) relative to the slopes crossing higher order replicas, until all slopes have been assigned a value ∆ max 1 (α), since by construction d 2 > d 1 , in the following, we only consider α > 1. The following procedure is obtained: Inside the set of all possible slopes (1, +∞), determine ∆ max 1 (α) for the subset of slopes α blocked by the first-order replica of the L(∆ 1 , ∆ 2 ) function, i.e., (k, h) = (0, 1). Notice that all the replicas (k, h) = (0, h) are blocked by replica (0,1) so they will be ignored in the following. II.
Determine the set of the non-blocked slopes after step I. III.
Determine the set of the non-blocked slopes after step III. V.
Determine the subset of slopes α blocked by next higher-order replica of the Determine the set of the non-blocked slopes after step V. VII.
Repeat steps V and VI until the set of non-blocked slopes is empty.
Conceptually, the presented procedure could be extended to N = 4 antenna elements, by operating in a 3D space, where the contours are replaced by surfaces, and the APP region resembles a Swiss cheese block. While the extension is straightforward, the analytical derivation would make the description uselessly cumbersome. Therefore, we avoid an explicit illustration of this case and leave it to the purview of interested reader.
As an example of the application of the N = 3 case, the proposed procedure is applied for the case with SNR = 20dB and P max = 10 −4 , resulting in L max = 8.1739. Since only slopes α > 1 are of interest, the APP region occupies only the portion above the bisector of the first and third quadrant of the plane (∆ 1 , ∆ 2 ), as shown in Figure 6, after remapping in the plane (∆ 1 , ∆ 2 − ∆ 1 ). As an example of the application of the = 3 case, the proposed procedure is applied for the case with = 20 and = 10 , resulting in = 8.1739 . Since only slopes > 1 are of interest, the APP region occupies only the portion above the bisector of the first and third quadrant of the plane (Δ , Δ ), as shown in Figure 6, after remapping in the plane (Δ , Δ − Δ ). Notice that the above procedure allowed us to determine the APP region in the plane (Δ , Δ ). By defining the extrema and of the angular sectors and , we can map the admissible region in the plane of the array parameters ( , ). The identified APP region will be the basis for the array design approach presented in the following section.

Three-Element NULA Design Strategy
In this section, we present the design strategy for the three-elements array aimed at Notice that the above procedure allowed us to determine the APP region in the plane (∆ 1 , ∆ 2 ). By defining the extrema u max and u int of the angular sectors U max and U int , we can map the admissible region in the plane of the array parameters (d 1 , d 2 ). The identified APP region will be the basis for the array design approach presented in the following section.

Three-Element NULA Design Strategy
In this section, we present the design strategy for the three-elements array aimed at low-cost low-weight sensors for DoA estimation. Based on the MSE approximation in Equation (7), it is clear that minimizing the CRB does not guarantee the best performance outside the asymptotic region, since for limited SNR, there might be a significant error contribution caused by the outliers that tend to concentrate around the DoAs, u m , where the beampattern exhibits its peaks. Moreover, Equation (7) shows that the values of P m act as scaling factors for the large error contributions (u m − u s ) 2 that add to the CRB.
Therefore, for a given SNR 0 , the only configurations that are of potential interest (i.e., provide small MSE) all have P m values below the maximum acceptable value, P max , and belong to the APP region. The procedure presented in Section 3 is therefore a useful tool to restrict the search for valuable array configurations and is used as an important part of the proposed design approach.
Inside the APP region, the optimum array configuration z = [0, d 1 , . . . , d N−1 ] is obtained, adopting a minimax criterion as the one that minimizes the maximum MSE over any desired range of source angles, U opt = −u opt , u opt , u opt = sin θ opt : U 0 can be either coincident with the whole sector of sources of interest U a (i.e., θ opt = θ int ), or it can be a smaller sector, where a higher accuracy is desired (i.e., θ opt < θ int ). This formulation is equivalent to: Since this array optimization approach is based on achieving the minimum MSE while keeping the probability of outliers under control, the obtained solution is referred to as the best outlier controlled array (BOCA). Its implementation diagram is sketched in Figure 7. The four blocks perform the following steps: (i) Identify the APP region in the plane (∆ 1 , ∆ 2 ), using the procedure in Section 3; (ii) Set the value of µ, the maximum value for u − u s , based on the sensor estimation characteristic and operational scenario; (iii) Map the APP region into the (d 1 , d 2 ) plane, based on the value of µ; (iv) Look for the array configuration (d 1 , d 2 ) providing the minimum of the MSE inside this region.
Notice that the APP region resulting from step (i) has been evaluated in the (∆ 1 , ∆ 2 ) plane. Recalling the definition of ∆ n in Equation (12), the APP region in Figure 6 can be remapped in the (d 1 , d 2 ) plane by scaling the maximum acceptable value ∆ max 1 (α) for every slope α by λ/µ max . Hence, the inter-element distances d n of the array configurations inside the APP region are given by (i) Identify the APP region in the plane (Δ , Δ ), using the procedure in Section 3; (ii) Set the value of , the maximum value for − , based on the sensor estimation characteristic and operational scenario; (iii) Map the APP region into the ( , ) plane, based on the value of ; (iv) Look for the array configuration ( , ) providing the minimum of the MSE inside this region. Notice that the APP region resulting from step (i) has been evaluated in the (Δ , Δ ) plane. Recalling the definition of Δ in Equation (12), the APP region in Figure 6 can be remapped in the ( , ) plane by scaling the maximum acceptable value Δ ( ) for The scaling in Equation (27) requires knowledge of µ max , which in turn requires defining the sectors U max , and U int . The bounds of these sectors largely depend on the considered application and purpose, as illustrated in the following section.
Equation (27) allows us to represent the APP region in the array inter-element distances axis system, i.e., in the (d 1 , d 2 ) plane, thus fulfilling step (iii). The final optimization algorithm (step (iv)) is based on a minimax criterion by choosing the array configuration that minimizes the maximum MSE obtained inside a given angular sector U opt = − sin θ opt , sin θ opt .
When compared to conventional, well-known array design strategies, such as the one presented in [28], the proposed approach offers three main advantages:

•
It allows for the relaxation of the half-wavelength quantization constraint on the inter-element distance. • It offers a high level of flexibility, as it allows us to distinguish the angular sector U max where the source can be located, the angular sector U int where the source is pursued, and the angular sector U opt , where the performance is optimized. • It allows us to obtain satisfactory DoA estimation performances, even when the SNR is in the threshold region, thanks to the constraint imposed on the maximum pairwise error probability.
In order to illustrate these advantages, in the next section, we assess the performances obtained by the BOCA in different case studies via numerical analysis.

Performance Assessment on Simulated Data
The purpose of this section is to test the proposed design approach against simulated data in order to prove its effectiveness. In addition, we compare it to the conventional design strategy presented in [28], where the array inter-element distances are subject to a λ/2 quantization. To make a fair comparison, we compare the proposed BOCA with a quantized array that satisfies the same constraint on the maximum pairwise probability. This can be achieved by selecting one of the quantized array configurations falling inside the APP region, as conducted for the BOCA. The selected array will be referred to as λ/2-best outlier controlled array (λ/2-BOCA). Due to the quantization constraint, we expect that the λ/2-BOCA will obtain an equal or higher MSE than the BOCA. Furthermore, there might be cases in which a λ/2-BOCA is impossible to locate, as all quantized configurations fall outside the APP region.
As mentioned in the previous sections, the APP region depends on both θ max and θ int , while the selected BOCA inside the APP region changes with θ opt . Since different values for θ max , θ int , and θ opt can be selected based on the requirements of different applications, we define two case studies that will be used for performance assessment: (i) Case A: θ max = 90 • , so that the source can be located at any possible DoA. Therefore, the possibility of having sources outside the angular sector U int of interest is considered. (ii) Case B: θ max = θ int , so that the possibility of having targets outside a given angular sector [−u max , u max ] is excluded. This can be useful when electromagnetic shielding is employed, or when the antennas are characterized by a very directive pattern.
Both case studies are considered in Sections 5 and 6, respectively dedicated to unconstrained and constrained design solutions. Notice that restricting θ max , as in case B, leads to smaller values of µ, and thus to larger APP regions, as shown in Equation (27). In turn, this allows us to choose longer array configurations, which exhibit better performances in terms of MSE.
To study the performance of the resulting BOCA, for both Case A and Case B, we evaluate the MSE of its DoA estimate inside the area of interest (−θ opt , θ opt ) as a function of the two remaining parameters, namely θ opt and θ int , assuming assigned values for SNR and P max . Since the MSE is variable in this region, a compact method to represent the performance is to average its value inside the sector U opt . This provides us with a single averaged MSE value for each choice of θ opt and θ int , since θ opt ≤ θ int , the angular sector U opt in which performances are studied, is a subset of the sector U int of the sources of interest, i.e.,: U opt ⊆ U int .
To show numerical results, we assume an SNR value of 20 dB, and we fix P max = 10 −4 . Figure 8a,b shows the map of the average MSE for case A and case B, respectively, for all the admissible values of θ int and θ opt . According to these results, the following observations apply:  To compare the performance of the BOCA and the /2-BOCA, it is convenient to define the ratio ( , ) between the two mean MSEs, namely: While the maximum (worst) MSE inside the sector for the BOCA is always To compare the performance of the BOCA and the λ/2-BOCA, it is convenient to define the ratio R θ opt , θ int between the two mean MSEs, namely: While the maximum (worst) MSE inside the sector U opt for the BOCA is always smaller compared to the λ/2-BOCA, the ratio R θ opt , θ int of the averaged MSEs over the sector U opt are not necessarily greater than the unity. This ratio allows us to quantify the average improvement achieved by the proposed approach over the conventional method, as a function of θ opt and θ int , over the whole area of interest. Figure 9a,b shows the R θ opt , θ int map for case study A and case study B, respectively. By observing these results, the following remarks apply: The maximum achievable ratio is ( , ) ≅ 3. This means that is at most about 3.5 times higher than ( , ) ).
(ii) As anticipated at the beginning of this section, there are cases in which no /2-BOCA can be found. For the selected parameters = 20 dB and = 10 , this occurs when choosing ≥ 80°. In these cases, the ratio shown in the figures has been conveniently saturated to its maximum value, namely 3.5, to denote that the improvement is not measurable. (iii) As previously mentioned, in case B, we obtain wider APP regions. Therefore, longer arrays are generally admissible, and the average MSE in case B is always smaller or than or equal to the average MSE in case A, for any given pair of and . (iv) Finally, while in case A, the best ratio values are obtained for smaller and , in case B, the best ratio values are obtained for = = ≈ 30°.
To complete the analysis, we show the detailed performance analysis for few specific parameter sets , , corresponding to the points of the , maps in Figure  9a,b, as shown with black circles, with the selected parameters , , and for each By observing these results, the following remarks apply: (i) The maximum achievable ratio is R θ opt , θ int ∼ = 3. This means that E MSE λ/2−BOCA θ opt , θ int is at most about 3.5 times higher than E MSE BOCA θ opt , θ int ). (ii) As anticipated at the beginning of this section, there are cases in which no λ/2-BOCA can be found. For the selected parameters SNR = 20 dB and P max = 10 −4 , this occurs when choosing θ a ≥ 80 • . In these cases, the ratio shown in the figures has been conveniently saturated to its maximum value, namely 3.5, to denote that the improvement is not measurable. (iii) As previously mentioned, in case B, we obtain wider APP regions. Therefore, longer arrays are generally admissible, and the average MSE in case B is always smaller or than or equal to the average MSE in case A, for any given pair of θ int and θ opt .
(iv) Finally, while in case A, the best ratio values are obtained for smaller θ int and θ opt , in case B, the best ratio values are obtained for θ max = θ int = θ opt ≈ 30 • .
To complete the analysis, we show the detailed performance analysis for few specific parameter sets θ opt , θ int , corresponding to the points of the R θ opt , θ int maps in Figure 9a,b, as shown with black circles, with the selected parameters θ max , θ int , and θ opt for each case study. In case study A, we focus on three subcases characterized by θ int = 30 • . However, while the average MSE of both the BOCA and the λ/2-BOCA monotonically increase with θ opt , their ratio shows a different behavior, namely R θ opt , θ int ≈ 1.25 at θ opt = 18 • and R θ opt , θ int ≈ 1.7 at both θ opt = 6 • and θ opt = 30 • . As the considered subcases are characterized by the same values of θ max , θ int , SNR, and P max , they all share the same admissible region, as represented in Figure 10. The location of the BOCA and /2-BOCA for the three subcases is rep a blue and a green circle, respectively, in the three subplots. With the sele rameters, the optimization procedure leads to the BOCA and /2-BOCA reported in Table 2.  The location of the BOCA and λ/2-BOCA for the three subcases is represented with a blue and a green circle, respectively, in the three subplots. With the selected set of parameters, the optimization procedure leads to the BOCA and λ/2-BOCA configurations reported in Table 2.    Figure 10) for | | ≤ 6°; however, for | | ≥ 12°, the BOCA is subject to the presence of a sidelobe, and its MSE shows a step increase, whereas the /2-BOCA is not subject to this effect. As clear from Figure 11b, for values of | | ≥ 12° , the BOCA configuration changes to = [0.00 0.66 5.22] , while the /2 -BOCA remains unchanged (see ○ markers in Figure 10). While the new BOCA still has a lower MSE than the /2-BOCA, its value is not as low as it was previously. When | | ≥ 22°, the array [0.00, 1.00, 4.50] is subject to the presence of a sidelobe, which provides a step increase in the MSE, so that the /2-BOCA configuration changes to / = [0.00, 0.50, 4.00], as apparent from Figure 11c. Its average MSE value increases, whereas the BOCA configuration is not subject to significant changes (see + markers in Figure 10). We also observe that when the MSE is too high, the simulated and theoretical MSE curves do not match exactly. This is because the union-bound approximation in Equations (7) and (8) is not tight.

Case Study A-Source
To complete the analysis, Figure 12 shows the maximum pairwise probability for each array configuration. The vertical dash-dot lines denote the area where the constraint must be guaranteed, while the horizontal dashed line represents the constraint value . Figure 12 illustrates that both the /2-BOCA and the BOCA satisfy the constraint on the maximum pairwise probability inside the angular sector , as we expected. However, as shown in Figure 11, the BOCA performs better inside the [− , ] angular sector in terms of MSE. The numerical MSE values were obtained through a Monte Carlo simulation, with N MC = 10 5 runs. The simulated array output x was generated according to the signal model in Equation (2), where the N random samples in the noise vector n were assumed to be complex valued and Gaussian distributed, namely: n ∼ CN 0, σ 2 n . The complexvalued amplitude A of the source baseband signal is a deterministic parameter, set so as to guarantee the required SNR condition, namely A = 2σ 2 n SNR. Figure 11a shows the MSE of z BOCA = [0.00 1. 30 5.84] to be lower than those of z λ/2−BOCA = [0.00, 1.00, 4.50] (see × markers in Figure 10) for |θ s | ≤ 6 • ; however, for |θ s | ≥ 12 • , the BOCA is subject to the presence of a sidelobe, and its MSE shows a step increase, whereas the λ/2-BOCA is not subject to this effect. As clear from Figure 11b, for values of |θ s | ≥ 12 • , the BOCA configuration changes to z BOCA = [0.00 0.66 5.22], while the λ/2-BOCA remains unchanged (see markers in Figure 10).
While the new BOCA still has a lower MSE than the λ/2-BOCA, its value is not as low as it was previously. When |θ s | ≥ 22 • , the array [0.00, 1.00, 4.50] is subject to the presence of a sidelobe, which provides a step increase in the MSE, so that the λ/2-BOCA configuration changes to z λ/2−BOCA = [0.00, 0.50, 4.00], as apparent from Figure 11c. Its average MSE value increases, whereas the BOCA configuration is not subject to significant changes (see + markers in Figure 10).
We also observe that when the MSE is too high, the simulated and theoretical MSE curves do not match exactly. This is because the union-bound approximation in Equations (7) and (8) is not tight.
To complete the analysis, Figure 12 shows the maximum pairwise probability for each array configuration. The vertical dash-dot lines denote the area where the constraint must be guaranteed, while the horizontal dashed line represents the constraint value P max . Figure 12 illustrates that both the λ/2-BOCA and the BOCA satisfy the constraint on the maximum pairwise probability inside the angular sector U int , as we expected. However, as shown in Figure 11, the BOCA performs better inside the −θ opt , θ opt angular sector in terms of MSE.  As opposed to case study A, in case B, the possibility that a source might be located outside the angular sector of interest is excluded, i.e., = . We focus on two subcases characterized by = 30°, as they seem to be the cases in which the best ratios , are obtained. As in case A, the different subcases are characterized by the same admissible region, as represented in Figure 13. With the selected set of parameters, the optimization procedure leads to the BOCA and /2-BOCA reported in Table 2.

Case Study B-Source Only Located within the Angular Sector
As opposed to case study A, in case B, the possibility that a source might be located outside the angular sector of interest is excluded, i.e., θ max = θ int . We focus on two subcases characterized by θ int = 30 • , as they seem to be the cases in which the best ratios R θ opt , θ int are obtained. As in case A, the different subcases are characterized by the same admissible region, as represented in Figure 13. With the selected set of parameters, the optimization procedure leads to the BOCA and λ/2-BOCA reported in Table 2.
Note that the decrease in θ max yields a larger admissible region. This allows us to choose longer array configurations, which were non-admissible in case A. Despite this, the λ/2-BOCA is still the one selected in the previous cases, confirming the reduced flexibility of the quantized approach. These considerations are further confirmed by Figure 14a,b and Figure 15a,b, where the MSE and the maximum pairwise error probability obtained using the different configurations is evaluated across a grid of θ values.
As a final consideration, we notice that the admissible regions shown in Figures 10 and 13, have very spiky shapes in certain areas. Selecting an array configuration inside one of those thin portions would mean requiring an installation accuracy that would be unlikely to be guaranteed in practice. Furthermore, these types of solutions are highly unstable, and lead to the variability of the ratio R θ opt , θ int with θ opt . Therefore, in the next section, we modify the optimization problem by adding additional constraints that might be required for the practical application of the proposed approach.

Interest ( = )
As opposed to case study A, in case B, the possibility that a source might be located outside the angular sector of interest is excluded, i.e., = . We focus on two subcases characterized by = 30°, as they seem to be the cases in which the best ratios , are obtained. As in case A, the different subcases are characterized by the same admissible region, as represented in Figure 13. With the selected set of parameters, the optimization procedure leads to the BOCA and /2-BOCA reported in Table 2. Note that the decrease in yields a larger admissible region. This allows us to choose longer array configurations, which were non-admissible in case A. Despite this, the /2-BOCA is still the one selected in the previous cases, confirming the reduced flexibility of the quantized approach. These considerations are further confirmed by Figure  14a,b and Figure 15a,b, where the MSE and the maximum pairwise error probability obtained using the different configurations is evaluated across a grid of values. Note that the decrease in yields a larger admissible region. This allows us to choose longer array configurations, which were non-admissible in case A. Despite this, the /2-BOCA is still the one selected in the previous cases, confirming the reduced flexibility of the quantized approach. These considerations are further confirmed by Figure  14a,b and Figure 15a,b, where the MSE and the maximum pairwise error probability obtained using the different configurations is evaluated across a grid of values. As a final consideration, we notice that the admissible regions shown in Figure 10 and Figure 13, have very spiky shapes in certain areas. Selecting an array configuration inside one of those thin portions would mean requiring an installation accuracy that would be unlikely to be guaranteed in practice. Furthermore, these types of solutions are highly unstable, and lead to the variability of the ratio , with . Therefore, in the next section, we modify the optimization problem by adding additional constraints that might be required for the practical application of the proposed approach. and Figure 13, have very spiky shapes in certain areas. Selecting an array configuration inside one of those thin portions would mean requiring an installation accuracy that would be unlikely to be guaranteed in practice. Furthermore, these types of solutions are highly unstable, and lead to the variability of the ratio , with . Therefore, in the next section, we modify the optimization problem by adding additional constraints that might be required for the practical application of the proposed approach.

Constrained Three-Element NULA Design Strategy
In this section, we consider two technological constraints that are likely to arise in real case scenarios, and we introduce a slight modification in the proposed strategy that will make it suitable for practical applications.
The first technological constraint is introduced based on the following observation: the APP regions shown in the previous sections always include solutions in which the inter-element distance is very small. In practical applications, depending on the signal wavelength, the physical size of the antennas, and the coupling effects that might occur, it might not be possible to position the antenna elements at the required proximity. Therefore, we consider the possibility to exclude the solutions characterized by inter-element distances shorter than a fixed value l from the admissible region.
The second technological constraint is related to installation accuracy. We define a finite installation accuracy δ, representing the maximum error tolerated in the installation.
Specifically, a given array z = 0 d 1 d 2 inside the APP region is considered admissible only if the four arrays 0 d 1 ± δ d 2 ± δ also lie within the admissible region. This excludes all the thinner areas from the admissible region.
The block diagram in Figure 16 illustrates the constrained optimization procedure. Clearly, for l = 0 and δ = 0, the constrained design strategy corresponds to the one in Equation (26).
To show the effect of the additional constraints on the results, let us consider a passive location system operating in the Wi-Fi band at f 0 = 2.447 GHz (λ = 0.1226 m). We assume that the employed commercial antennas have a size l ≈ 0.13 m = 1.1λ, constraining the minimum inter-element distance. Two different values are considered for δ, namely δ = 0.005λ = 6.1·10 −4 m (micrometric positioning) and δ = 0.05λ = 6.1·10 −3 m (manual positioning). Figure 17a,b shows the admissible region obtained by including the technological constraints for the two cases, assuming SNR = 20 dB, P max = 10 −4 and θ max = θ int = 30 • . The blue-colored areas are the portions of the APP region excluded by the constraint on the minimum inter-element distance, while the red-colored areas are the portions of the APP excluded to consider installation accuracy tolerance. As visible from Figure 17, the installation accuracy constraint removes thinner areas from the APP region. These areas contain array configurations with highly unstable performance, for which a slight modification in the inter-element distance causes the array configuration to fall outside the APP region.
it might not be possible to position the antenna elements at the required proximity. Therefore, we consider the possibility to exclude the solutions characterized by inter-element distances shorter than a fixed value from the admissible region.
The second technological constraint is related to installation accuracy. We define a finite installation accuracy , representing the maximum error tolerated in the installation. Specifically, a given array = 0 ̅ ̅ inside the APP region is considered admissible only if the four arrays 0 ̅ ± ̅ ± also lie within the admissible region. This excludes all the thinner areas from the admissible region.
The block diagram in Figure 16 illustrates the constrained optimization procedure. Clearly, for = 0 and = 0 , the constrained design strategy corresponds to the one in Equation (26).  Figure 17a,b shows the admissible region obtained by including the technological constraints for the two cases, assuming = 20 dB, = 10 and = = 30°. The blue-colored areas are the portions of the APP region excluded by the constraint on the minimum inter-element distance, while the red-colored areas are the portions of the APP excluded to consider installation accuracy tolerance. As visible from Figure 17, the installation accuracy constraint removes thinner areas from the APP region. These areas contain array configurations with highly unstable performance, for which a slight modification in the inter-element distance causes the array configuration to fall outside the APP region.
(a) (b) As visible in Figure 17, the technological constraints lead to significant changes in the APP region, so that most of the array configurations considered in case studies A1-A3, and B1 and B2, are no longer admissible. Therefore, for each choice of , , and , new BOCA and /2-BOCA configurations must be identified. For example, considering = 1.1 and = 0.05 , as shown in Figure 17b, and assuming = = = 30° , As visible in Figure 17, the technological constraints lead to significant changes in the APP region, so that most of the array configurations considered in case studies A1-A3, and B1 and B2, are no longer admissible. Therefore, for each choice of θ max , θ int , and θ opt , new BOCA and λ/2-BOCA configurations must be identified. For example, considering l = 1.1λ and δ = 0.05λ, as shown in Figure 17b, and assuming θ max = θ int = θ opt = 30 • , the new BOCA configurations become z BOCA = [0 1.96 6.82]λ and z λ/2−BOCA = [0 1.5 3.5]λ. We also notice that the obtained BOCA configuration z BOCA is characterized by a very low value of min S pq = 4.74, as evaluated following the procedure in [31], and it would have not been considered admissible based on the maximum error criterion. However, z BOCA is guaranteed to satisfy the constraint P m < P max , since it belongs to the APP region. Therefore, in this case, the min S pq value is not a good proxy for robustness to estimation ambiguities, since our ML-based admissibility criterion allows us to find a robust array configuration with low min S pq . This is further confirmed by observing the ratio R θ opt , θ int between the MSEs of the z BOCA and that of the z λ/2−BOCA . Since we have R θ opt , θ int = 3.1510, the BOCA design solution ultimately achieves an MSE which is about one-third of the MSE obtained by the λ/2-BOCA. In conclusion, z BOCA is indeed a good alternative to both the λ/2-BOCA and to the BOCA designs obtained in the case studies without technological constraints.
Depending on the value of the constraints and on the extrema of the considered angular sectors, the modified APP region might become so restrictive that no λ/2-BOCA configurations could be found. In such cases, the ratio R θ opt , θ int is saturated at its maximum value, namely R θ opt , θ int = 3.5. Figure 18a,b shows the contour plots of R θ opt , θ int , obtained considering case studies A and B, respectively, and assuming an installation accuracy tolerance of δ = 0.005λ. As is apparent, the additional design constraints degrade the estimation accuracies of both the λ/2-BOCA and the BOCA configurations. Therefore, the ratio R between the two MSEs still assumes values between R ∼ = 1 and R ∼ = 3.5. installation accuracy tolerance of = 0.005 . As is apparent, the additional design constraints degrade the estimation accuracies of both the /2-BOCA and the BOCA configurations. Therefore, the ratio between the two MSEs still assumes values between ≅ 1 and ≅ 3.5. By increasing the installation accuracy tolerance to = 0.05 , we obtain the contour plots of , , as shown in Figure 19a,b, relative to cases A and B, respectively. By increasing the installation accuracy tolerance to δ = 0.05λ, we obtain the contour plots of R θ opt , θ int , as shown in Figure 19a,b, relative to cases A and B, respectively. By increasing the installation accuracy tolerance to = 0.05 , we obtain the contour plots of , , as shown in Figure 19a,b, relative to cases A and B, respectively.
(a) (b) Figure 19. The 2D improvement maps for the constrained case with = 0.05 : (a) case study A; (b) case study B.
As visible from these contour plots, the effect of the higher tolerance can provide a significant change in the array design. Specifically, most of the , plane is saturated to 3.5. Based on this observation, we conclude that the proposed approach is always able to provide a valid array configuration, while there is no /2-BOCA solution that satisfies all the imposed constraints at the same time. As visible from these contour plots, the effect of the higher tolerance can provide a significant change in the array design. Specifically, most of the θ opt , θ int plane is saturated to 3.5. Based on this observation, we conclude that the proposed approach is always able to provide a valid array configuration, while there is no λ/2-BOCA solution that satisfies all the imposed constraints at the same time.
In conclusion, the results obtained in case study B, assuming δ = 0.005λ, are also used to validate the obtained design strategy in a real-case scenario using experimental data, as reported in the following section.

Experimental Results
To further support the effectiveness of the design procedure proposed in the previous sections, an experimental trial was conducted in a parking area.
For this experiment, we assumed θ max = θ int = 45 • . The design of the array has been carried out to guarantee P m < P max = 10 −4 when the operative SNR value is set to SNR 0 = 20 dB. Finally, the minimum inter-element distance has been constrained to be l = 1.1λ, as in the previous section, with δ = 0.005λ. With these parameters, we obtain the admissible region in Figure 20.
To carry out the experimental acquisition, we exploited a four-channel National Instruments NI USRP-2955 board, operating in the 2.5 GHz Wi-Fi frequency band. To collect samples simultaneously with the two arrays, we must operate with two configurations that have two elements in common. For this purpose, we selected a non-optimal array, addressed in the following as Almost-BOCA array z ABOCA = [0 2 5.35]λ, that has one interelement distance in common with the λ/2-BOCA, z λ/2−BOCA = [0 2 4.5]. This allowed us to make the acquisitions with a single 4-elements array with z = [0 2 4.5 5.35], with the individual element connected channels 0,1,2,3 of the 4-channel USRP acting as the coherent receiver. With this arrangement, the snapshots of the λ/2-BOCA array were obtained by using only the samples received by channels 0,1, and 2, whereas those of the Almost-BOCA array were obtained by using only the samples received by channels 0,1, and 3. Moreover, in this way, the snapshots of the two arrays are coherent and simultaneous. Note that the resulting array z ABOCA provides only slightly degraded performance with respect to the BOCA.
For this experiment, we assumed = = 45°. The design of the arr carried out to guarantee < = 10 when the operative value is se 20 dB. Finally, the minimum inter-element distance has been constrained to as in the previous section, with = 0.005 . With these parameters, we obtain sible region in Figure 20. To carry out the experimental acquisition, we exploited a four-channel N struments NI USRP-2955 board, operating in the 2.5 GHz Wi-Fi frequency ban samples simultaneously with the two arrays, we must operate with two con that have two elements in common. For this purpose, we selected a non-op addressed in the following as Almost-BOCA array individual element connected channels 0,1,2,3 of the 4-channel USRP acting a ent receiver. With this arrangement, the snapshots of the /2-BOCA array we by using only the samples received by channels 0,1, and 2, whereas those of BOCA array were obtained by using only the samples received by channels Moreover, in this way, the snapshots of the two arrays are coherent and sim Note that the resulting array provides only slightly degraded perfor respect to the BOCA.  Figure 21a shows a sketch of the acquisition geometry. The green markers represent the element of the λ/2-BOCA, while the blue markers indicate the elements of the Almost-BOCA. The two sensors in the middle are common between the two arrays. The fourelement array that includes both the Almost-BOCA and the λ/2-BOCA solutions was realized by deploying four Wi-Fi antennas on a plastic support mounted on a tripod. Figure 21b shows this receiving system set in the origin of the reference system. The Wi-Fi access point (AP) in Figure 21c is used as the transmitting source, which has been initially placed at 20 m with DoA θ s = 0 • from the receiving antennas and then moved to θ s = 30 • with respect to the boresight of the array. This was set to transmit the Wi-Fi beacons with a beacon interval of 3ms at f 0 = 2.447 GHz, (namely, using wavelength λ = 0.1226 m).
The expected performance for SNR = 20 dB in terms of both MSE and P out is reported in Figure 22a,b, as obtained from the theoretical framework in Equations (7)-(9) for both the Almost-BOCA and the λ/2-BOCA. The following remarks apply: When the source is located at θ s = 0 • , the MSE obtained with the λ/2-BOCA should be about 2 times higher than the MSE obtained with the Almost-BOCA, while the P out is about 4 times lower.
When the source is located at θ s = 30 • , the MSE obtained with the λ/2-BOCA should be about 1.5 times higher than the MSE obtained with the Almost-BOCA, while the P out is about 3 times lower. element array that includes both the Almost-BOCA and the λ/2-BOCA solutions was realized by deploying four Wi-Fi antennas on a plastic support mounted on a tripod. Figure  21b shows this receiving system set in the origin of the reference system. The Wi-Fi access point (AP) in Figure 21c is used as the transmitting source, which has been initially placed at 20 m with DoA = 0° from the receiving antennas and then moved to = 30° with respect to the boresight of the array. This was set to transmit the Wi-Fi beacons with a beacon interval of 3ms at = 2.447 GHz, (namely, using wavelength = 0.1226 m). The expected performance for = 20 dB in terms of both MSE and is reported in Figure 22a,b, as obtained from the theoretical framework in Equations (7)- (9) for both the Almost-BOCA and the λ/2-BOCA. The following remarks apply: When the source is located at = 0°, the MSE obtained with the /2-BOCA should be about 2 times higher than the MSE obtained with the Almost-BOCA, while the is about 4 times lower.
When the source is located at = 30°, the MSE obtained with the /2-BOCA should be about 1.5 times higher than the MSE obtained with the Almost-BOCA, while the is about 3 times lower.
The experimental results have been obtained by simultaneously collecting long sequences of samples with the 4 elements when the AP is active. These data are received with a measured = 23 dB. White Gaussian noise was added to the data to emulate lower values, specifically between 14 dB and 23 dB, with a 0.5 dB step. By averaging the received samples, the resulting MSE and curves are obtained and reported in Figure 23a,b, respectively, as a function of the value and for both considered DoAs. By observing Figure 23, we note that:

•
For both the considered DoAs, the Almost-BOCA solution outperforms the λ/2 -BOCA in the asymptotic region, since a longer array is used; • For = 0°, the Almost-BOCA solution outperforms the λ/2-BOCA, both in terms of MSE and for ≥ 16 dB; • For = 30°, the improvement is lower as predicted by theory, and starts at a higher ( ≥ 19 dB). The experimental results have been obtained by simultaneously collecting long sequences of samples with the 4 elements when the AP is active. These data are received with a measured SNR = 23 dB. White Gaussian noise was added to the data to emulate lower SNR values, specifically between 14 dB and 23 dB, with a 0.5 dB step.
By averaging the received samples, the resulting MSE and P out curves are obtained and reported in Figure 23a,b, respectively, as a function of the SNR value and for both considered DoAs. By observing Figure 23, we note that:

Conclusions
DoA estimation of narrowband signals is a ubiquitous task in several civil surveillance applications, such as passive coherent location systems, passive sonar, or passive radar. The main challenge in these applications is to comply with the typical low-cost requirements involving not only the economic cost, but also the computational load. Additionally, the DoA estimation accuracy is usually degraded by the presence of outliers, which occur as a result of the poor conditions characterizing low-cost systems. To achieve the required features of limited complexity and system lightness, the number of antenna elements must be kept sufficiently low. With this perspective, NULA configurations represent a viable solution to reduce the number of receiving elements without compromising performance.
The design strategy presented in this paper is set within the context of low-cost sensors applications, as it allows us to identify the so-called best outlier-controlled NULA configuration, achieving satisfactory DoA estimation performance by exploiting just three antenna elements. As a matter of fact, the BOCA solution allows for minimizing the MSE, while keeping the outlier probability under reasonable control.
The design procedure is based on the following two steps. First, an admissible pairwise probability region is identified, exploiting the algorithm presented in Section 3. The APP region contains all the arrays for which the outlier probability is kept under control. Second, the BOCA configuration achieving the minimum MSE is located inside the APP region. It has also been shown that the APP region can be easily modified to incorporate practical design constraints, such as the minimum inter-element distance and the installation accuracy tolerance.
The proposed design strategy has been compared to the one presented in [28], which represents our benchmark. Therein, the array inter-element distances are quantized to half-wavelength. To carry out a fair analysis, we compared the estimation accuracy achieved by the BOCA with the one achieved by the best outlier-controlled array with the half-wavelength quantization constraint, namely the /2-BOCA. As shown via both theoretical and numerical analyses, the BOCA solutions generally achieve improved estimation accuracies compared to the /2-BOCA estimations, with an improvement depending on how the angular sectors , , and are chosen. The effectiveness of the proposed design procedure has also been verified using experimental data, confirming that the optimized array outperforms the /2 -quantized benchmark array. The comparison with a recent three-element NULA design approach to control ambiguity outliers is also

•
For both the considered DoAs, the Almost-BOCA solution outperforms the λ/2-BOCA in the asymptotic region, since a longer array is used; • For θ s = 0 • , the Almost-BOCA solution outperforms the λ/2-BOCA, both in terms of MSE and P out for SNR ≥ 16 dB; • For θ s = 30 • , the improvement is lower as predicted by theory, and starts at a higher SNR (SNR ≥ 19 dB).
As apparent, the experimental results show that the non-λ/2-spaced array always outperforms its λ/2 counterpart for the design value SNR 0 = 20 dB and for higher SNR values. This confirms that the design procedure effectively allows us to design an array that provides a better performance than the λ/2-BOCA in a realistic experimental scenario.
We recall that a slightly different and better performing BOCA solution would be selected if we did not force it to have one inter-element distance in common with the λ/2-BOCA. Therefore, higher improvements are expected if the design strategy is applied with its entire set of degrees of freedom.

Conclusions
DoA estimation of narrowband signals is a ubiquitous task in several civil surveillance applications, such as passive coherent location systems, passive sonar, or passive radar. The main challenge in these applications is to comply with the typical low-cost requirements involving not only the economic cost, but also the computational load. Additionally, the DoA estimation accuracy is usually degraded by the presence of outliers, which occur as a result of the poor SNR conditions characterizing low-cost systems. To achieve the required features of limited complexity and system lightness, the number of antenna elements must be kept sufficiently low. With this perspective, NULA configurations represent a viable solution to reduce the number of receiving elements without compromising performance.
The design strategy presented in this paper is set within the context of low-cost sensors applications, as it allows us to identify the so-called best outlier-controlled NULA configuration, achieving satisfactory DoA estimation performance by exploiting just three antenna elements. As a matter of fact, the BOCA solution allows for minimizing the MSE, while keeping the outlier probability under reasonable control.
The design procedure is based on the following two steps. First, an admissible pairwise probability region is identified, exploiting the algorithm presented in Section 3. The APP region contains all the arrays for which the outlier probability is kept under control. Second, the BOCA configuration achieving the minimum MSE is located inside the APP region. It has also been shown that the APP region can be easily modified to incorporate practical design constraints, such as the minimum inter-element distance and the installation accuracy tolerance.
The proposed design strategy has been compared to the one presented in [28], which represents our benchmark. Therein, the array inter-element distances are quantized to halfwavelength. To carry out a fair analysis, we compared the estimation accuracy achieved by the BOCA with the one achieved by the best outlier-controlled array with the halfwavelength quantization constraint, namely the λ/2-BOCA. As shown via both theoretical and numerical analyses, the BOCA solutions generally achieve improved estimation accuracies compared to the λ/2-BOCA estimations, with an improvement depending on how the angular sectors U max , U int , and U opt are chosen. The effectiveness of the proposed design procedure has also been verified using experimental data, confirming that the optimized array outperforms the λ/2-quantized benchmark array. The comparison with a recent three-element NULA design approach to control ambiguity outliers is also included, which emphasizes the capability of our proposed model to accurately characterize and control the probability of ambiguities.
As a future research scope, it might be interesting to extend the design strategy to four-element and five-element NULA configurations in order to improve the estimation accuracy without significantly increasing the number of receiving channels. In such cases, alternative design strategies, such as nested or coprime array configurations, may also provide a useful benchmark. Furthermore, non-linear array configurations, such as circular arrays, could be investigated, as they would allow us to estimate the DoA regarding both azimuth and elevation. Lastly, the DoA estimation problem could be extended to multisource scenarios, which would be useful not only when multiple sources of interest are present, but also to take any multipath effect into proper account.