Non-Uniform Warping Sampling for Data Reduction in Planar Array Diagnostics

The problem of detecting defective turned-off elements in antenna arrays from near-field measurements is addressed. In particular, the focus here is to reduce the number of measurements in order to positively affect the acquisition time. Such an issue is achieved by adopting the recently developed warping sampling method. Two commonly antenna diagnostics methods, i.e, the Back Transformation Method (BTM) and the Matrix Method (MM), are considered in view of this new sampling strategy and compared to the usual half-wavelength sampling. In particular, in order to identify the fault locations, outcomes returned by BTM and MM undergo a detection step based on a cell-averaging CFAR (CA)-CFAR technique borrowed from the radar literature. It is shown that the warping sampling method provides performance close to the uniform half-wavelength one with a reduced number of data. Numerical simulations are carried out in order to verify the results with different fault layouts and tapered currents.


I. INTRODUCTION
Array antennas are widely used in many applications such as radar, automotive, wireless communication etc. They consist of a number of radiating elements which, depending on the applications, can sometimes be very large.
An essential step in array development is the diagnostic stage, which aims at checking if the array complies with the design specifications. This entails, for example, looking for defective (faulty) elements, which from a mathematical point of view consists in solving an inverse source problem, i.e., estimating the excitation coefficients from field measurements [1]- [7].
For planar arrays diagnostics, the most commonly used method is by far the Back Transformation Method (BTM) [8]. The Matrix Method (MM) [9] is also used often as it has greater flexibility and allows to deal with generally shaped array and measurement apertures.
The associate editor coordinating the review of this manuscript and approving it for publication was Qi Luo .
Classically, measurements are collected over an aperture larger than the array support (to control truncation error [10]) over a uniform λ/2 grid. This naturally matches the standard FFT routines used to implement the BTM. However, this can yield a lot of data resulting in prolonged measurement collection time, which can be a serious inconvenience. Hence, reducing the required data would speed up the diagnostics stage.
To this end, approaches based on compressed sensing (CS) [11]- [15] have been recently explored. Indeed, by formulating the problem as the search for defective elements and by furthermore assuming that their number is low, CS allows to retrieve a sparse solution by using a number of data below the Nyquist rate [16], [17] by running a l 1 minimization. In this framework, ensuring the RIP condition is crucial. Sensing matrices, which statistically verify the RIP condition, can be built by using random matrix theory. This can be done, for example, by randomly selecting a subset of sampling positions. This method is known as random sampling (RS). However, the limitation of such an approach is that the RIP of the sensing matrix is assured in probability. Without a proper rule to select the elements of the available ensemble, there exists a probability, even if small, that the inversion might be unstable, causing erroneous fault detection. In [15], a method to pick up a specific instance from the ensemble of possible sampling positions assuring the RIP property is proposed. However, it allows to obtain a deterministic sampling (DS) strategy for far-field measurement. Results referring to near-field measurements are not available. Also, to properly set the number of data, CS relies on the knowledge of the unknown sparsity.
The number of degrees of freedom of the radiated field (NDF) [18], [19], depending on the size of the source and the measurement aperture, can be much lower than the number of points returned by the λ/2 sampling [20]. This fact has been indeed exploited in [21], where the so-called warping method was introduced and used to derive a new deterministic near-field sampling strategy. This allowed for a remarkable reduction of data points by keeping accuracy close to the λ/2 sampling while estimating the radiation pattern [22]. This result can be in principle useful for CS as well, since the points it returns actually represent the baseline from which further data reduction can be achieved if unknown sparsity information is available. On the other hand, it allows for directly exploiting MM or BTM. Indeed, since the data points result non-uniformly arranged, BTM requires a preliminary interpolation stage. Accordingly, despite a great data reduction (as compared to the λ/2 sampling) standard inversion methods can be employed to achieve diagnostics. This is interesting since we do not need to care about RIP or to run some optimization to achieve inversion. Also, noise propagation can be better controlled because of the involved standard inversion procedures. Finally, no a priori information about the unknown sparsity (equivalent to knowing the fault percentage) is required.
Accordingly, in this contribution, we employ BTM and MM for achieving array diagnostics, by considering as faulty elements the ones that are turned off. BTM and MM return the elements that actually populate the array (i.e., the ones that are working properly), hence the defective ones can be identified as the nulls (actually minima) of the reconstructions. This can be difficult to achieve because of the noise and the filtering introduced by the inversion procedures. To cope with this limitation, faulty elements are emphasized by subtracting from the reconstructed coefficients the actual ones (assumed known) projected on the same subspace the reconstructions belong to. Also, a fault detection step follows. It consists in introducing an adaptive detection threshold determined according to the CFAR technique borrowed from radar literature [23].
Summarizing, this contribution is based on two main ingredients: the warping sampling strategy developed in [22], which allows to perform diagnostics by a significant reduction of data points as compared to the λ/2 sampling, and the application of a detection strategy to highlight the turned-off elements in the array. An extensive numerical analysis is included in order to assess the achievable performance.

II. ARRAY DIAGNOSTICS FORMULATION
Consider the antenna under test (AUT) sketched in Fig. 1.
of the x − y plane, according to a uniform rectangular grid. Denote as r n = (x n , y n , 0) and c n the element positions and their excitation coefficients respectively and as f the element factor.
The field radiated by such an antenna is collected over a planar observation domain r 0 ∈ OD = [−X 0 , X 0 ]×[−Y 0 , Y 0 ] located at z = z 0 > λ, with λ being the wavelength, so that only propagating waves are relevant.
This radiation problem can be described by the following semi-discrete radiation operator R(r 0 , r n ) = |r 0 − r n |, θ(r 0 , r n ) and φ(r 0 , r n ) the relative polar angles between the measurement point r 0 and the n-th element position r n , and k is the wavenumber. Measurements are collected by a probe antenna, whose plane-wave spectrum is generally flat enough to be ignored. Otherwise, a probe compensation procedure can be applied. Here, we straightway assume the ideal probe with a Dirac pulse response. Moreover, the probe is assumed linearly polarized and collecting only tangent field components E t (r 0 ) = E(r 0 ) ·t, witht =x ort =ŷ.
The array diagnostics problem aims to reconstruct the excitation coefficient vector c from field measurements. The back transformation method and the matrix method are two commonly used approaches to pursue such an objective and are briefly described below. VOLUME 10, 2022

A. BACK TRANSFORMATION METHOD
The back transformation method is a diagnostic technique based on field plane-wave spectrum (PSW) expansion.
SayÊ t (k x , k y , z 0 ) the field PSW at the measurement aperture z = z 0 , with k x and k y being the spatial spectral variables. Then, after compensating the propagation term and the element factor, an estimate of the excitation coefficients c n is obtained aŝ with f t = f ·t being the tangent component of the element factor. It is seen that estimating the excitation coefficients via BTM entails computing a couple of Fourier transformations, one for field PWS evaluation and one for achieving backpropagation. This can be achieved by employing an FFT routine, which makes the method appealing because of its quick implementation. Usually, the field is sampled at λ/2, λ being the wavelength, which is equivalent to restricting the band to the so-called visible domain, i.e., evanescent contribution is filtered out. Note that restricting the spatial spectrum implicitly regularizes the inverse problem of retrieving c n from field measurements E t (k x , k y , z 0 ), which is indeed illposed. This is easily seen by looking at the exponential term e jk z z 0 , which diverges for k 2 x + k 2 y > k 2 . However, (5) does not return the least square solution because of the finite size of the measurement aperture. This of course negatively affects the performance achievable in the reconstructions (truncation error) [10].

B. MATRIX METHOD (MM)
Since field data are collected over a finite discrete set of points belonging to OD, instead of the operator in (2), what one should actually consider is the matrix operator where E t is the numerical vector consisting of the samples of E t collected over the M measurement points r 0m and A is a matrix whose entries are given by where R mn = |r 0m − r n | and f tmn is the tangent component of the element factor linking the element at r n to the observation point at r 0m . The MM, hence, entails recovering the excitation vector c by solving the matrix equation (6), typically in the least square sense since M > N and because of noise and uncertainty. Clearly, MM is more computationally demanding than BTM, with the actual computational cost being related to the employed solution algorithm. Indeed, direct or iterative inversion schemes can be exploited. However, MM is more flexible since it allows to just as easily deal with measurement aperture and array antenna of general shapes (i.e., not necessarily planar). Data collected over a non-uniform grid can be addressed as well.
Regardless of the inversion method one may want to use, said {σ n , u n , v n } N n=0 the singular system of A, the retrieved excitation vector can be formally written aŝ where < ·, · > denotes the scalar product, N is the noise vector that corrupts data and W n is a filtering window that depends on the adopted inversion scheme and regularizes the problem. Regularization depends on the noise level and, of course, impacts on the achievable performance. However, in general, MM can allow a for better performance than BTM, since the latter ''imposes'' regularization regardless of the noise level.
Related to the previous question there is the strategy adopted to design the measurement set-up in terms of the number and the location of the measurement points. This crucial aspect is where we mainly contribute in this paper as detailed in the next section.

III. FIELD SAMPLING
As mentioned above, retrieving the excitation coefficients entails solving a finite dimensional ill-posed linear inverse problem, which basically inherits ill-posedness from the related continuous inverse source problem of which it is the discrete counterpart. More in detail, it has been shown that the ''continuous'' radiation operator (i.e., operator in (1) acting on L 2 (SD)) presents a distinctive singular value behaviour. Indeed, the singular values are rather constant till a critical index beyond which they experience an exponentially fast decay. This critical number is the so-called number of degrees of freedom (NDF) of the problem and basically means that the dimension of the radiated field space (i.e., the data space) is essentially of finite dimension. Accordingly, the number of spatial measurements, and their positions, should be dictated by the NDF (indeed slightly greater) so as to approximate such a finite dimensional space.
By the usual uniform λ/2 sampling, M > N > NDF. This is because in order to reduce truncation error the measurement aperture is commonly set larger than the array size. Also, note that typically NDF < N , that is lower than the number of elemental radiators uniformly deployed at a λ/2 step over the array planar aperture. It has been shown that the NDF → N only when the measurement aperture becomes unbounded, which is obviously unfeasible. More in detail, the NDF was analytically estimated in [20] as with [·] being the operator that retains the least greatest integer. From (9) it is easily seen that when X 0 , Y 0 → ∞ then NDF → 16X S Y S , which is exactly N , that is the number of radiators populating the array. Therefore, when M > N > NDF, the singular values of the discrete operator (6) run beyond the NDF and hence very low singular values may lead to instability in the reconstruction (ill-posedness) [24]. This is circumvented in the BTM, since it implicitly regularizes the problem, whereas in the MM one has to properly choose the regularizing sequence W n . Since the radiated field space has essentially dimension NDF, it is natural to try to exploit only a number of measurements of the same order as the NDF. However, the measurement points should be set so that the singular values of the discrete operator A approximate the ones of G in the region preceding the abrupt decay. If this is possible, a great data reduction is achieved with consequent data collection time saving, especially for large or very large (in terms of wavelength) array antennas.
This crucial point has been recently successfully pursued in [21], [22]. There, a field sampling strategy, that requires a number of measurements M only slightly greater than the NDF reported in (9), was derived. This approach is based on a proper reformulation of the radiation operator that highlights the radiated field as a spatially varying bandlimited function, that is, as a function whose bandwidth depends on the observation point r 0 . More in detail, it is shown that upon introducing suitable transformations that ''warp'' the original observation variable r 0 , the radiated field can be expressed as a classical band-limited function and thus the Shannon sampling theorem can be employed. The transformations are shown in eq. (10), whereas the theoretical details are reported in [22].
where α x and α y are given by Indeed, ν is an oversampling factor so that M = ν 2 NDF and the sin terms allow to achieve a spatially varying oversampling so as to deploy more points where required (see [22]). Accordingly, the following sampling series can be exploited to represent the field Finally, the sampling points are evaluated by solving with s, l being integer numbers. It is seen that, since (ξ x , ξ y ) are non-linearly linked to (x 0 , y 0 ), uniform sampling in (ξ x , ξ y ) becomes non-uniform in r 0 . This could be expected in view of the spatially varying behaviour of the radiation operator in the near-field. Nonetheless, as remarked above, a considerable reduction of data points is achieved. Fig. 2 shows the ratio between the number of measurements M = ν 2 NDF returned by non-uniform sampling and M λ/2 , the number of measurements point corresponding to the uniform λ/2 sampling, as function of the observation domain size X 0 = Y 0 for different values of X S = Y S and by choosing ν = 1.2. As can be seen, such a ratio is always lesser than 1. Hence, the proposed sampling strategy allows to decrease the measurement number compared to the classical λ/2 sampling. In particular, the data reduction (with respect to M λ/2 ) becomes dramatic as the array size (and consequently the measurement aperture) increases.
Clearly reducing the number of measurement points is positive since data acquisition time is reduced as well. Also, MM enjoys a lighter computational burden since it is required VOLUME 10, 2022 to deal with a smaller matrix. BTM, instead cannot be employed directly. However, the series in (12) allows to interpolate the field over a λ/2 grid and then BTM can still employed.
The crucial question is to see how data reduction, allowed by the non-uniform sampling, affects performance in the reconstructions, compared to the λ/2 sampling. To this end, in the next section some examples of array diagnostics are provided.

IV. NUMERICAL EXAMPLES
We consider an array of N = 2809 dipoles directed along the x axis, uniformly spaced at λ/2 and arranged over the array support X s = Y s = 13λ. The x component of the radiated field is collected over a planar observation domain with X 0 = Y 0 = 25λ and z 0 = 8λ. The field is generated by the Phased-array toolbox of MATLAB. In particular, since such a toolbox accepts voltage taper, the desired excitation currents have been scaled by the antenna input impedance provided by the MATLAB toolbox itself.
Two different excitation distributions are considered: uniform, with c n = 1 ∀n, and ann = 4 Taylor distribution with side-lobe level SLL = 20dB. For comparison purposes, data are collected by exploiting both the usual λ/2 sampling and the non-uniform one returned by eq. (14). The first sampling scheme requires M λ/2 = 10201, whereas the non-uniform one M = 3600 (i.e., only 35% of M λ/2 ).
We focus on detecting defective elements. To this end, a fraction of the overall excitation coefficients is turned off according to a uniform probabilistic law. Both the MM and BTM are employed to reconstruct the excitation coefficients and each of them is checked for both the sampling schemes under comparison. Therefore, we have actually four reconstruction schemes, with BTM proceeded by interpolation (as mentioned above) when applied to the data obtained over the non-uniform grid. Finally, data are corrupted by a zero mean complex white Gaussian noise whose signal to noise ratio (SNR) is Fig. (3) shows two random faulty deployments for the case of uniform and Taylor excitation coefficients. In particular, the number of faulty elements was set at 3%N (N fault = 87) and the field data corrupted by noise with SNR = 20dB. The corresponding reconstructed excitation coefficients are shown in Fig.4 and Fig. 5. In particular, the MM reconstructions have been obtained by employing a truncated singular value decomposition (TSVD) inversion scheme and retaining the singular values above −20dB the maximum one. In both figures the top rows refer to MM and the bottom ones to BTM, instead the first columns have been obtained by considering the uniform λ/2 sampling whereas the right column using the non-uniform sampling. It can be appreciated that, from the considered singular value truncation threshold, the MM and the BTM reconstructions look very similar. What is more, this holds true for both the sampling schemes, in spite of the reduction in the number of measurements of the non-uniform sampling. This is sufficient to state that the non-uniform sampling scheme actually captures the same information as the λ/2 sampling. Also, it can be noted that most of the faulty elements are already clearly detectable, especially for the uniform excitation. However, for the Taylor's taper, discerning from faults and excitation coefficients becomes increasingly hard in the region where the latter are very low. This issue can be remedied by employing a suitable detection strategy.

V. FAULT DETECTION
To begin with, first the difference excitation coefficients are built. In particular, the difference coefficients are defined as c = Pc −ĉ, whereĉ are the reconstructed coefficients and Pc are the actual excitation coefficients (assumed known) that have undergone the same filtering (P represents the corresponding projection operator) as due to the reconstruction procedure. For example, for the MM (similar arguments apply for BTM), it yields (16) where N T is the TSVD truncation index and I fault is the set of positions where faults occur. It is noted that knowing Pc is equivalent to knowing the field radiated by the array when all elements work correctly, as commonly done in many papers concerning array diagnostics [11]- [15]. Also, eq. (16) can be expressed in terms of the point-spread function (psf) of the adopted inversion scheme as which is useful for later discussion. Detection is then achieved by turning c into a binary vector populated by 1 where | c(n)| 2 ≥ A th , A th being the detection threshold, and 0 elsewhere.
The choice of the detection threshold is a crucial issue. Indeed, it must be set in order to establish a proper balance between the probability of false alarm (PFA) and the probability of detection (PD). In the case of a single fault, and under the assumed working condition, analytical results that link A th and PFA and PD can be established, as done in [25]. However, when multiple faults occur, because of the linkage between the reconstructions of the different faults, those results can be expected to hold only when the faults are well separated. The point is that one does know which is the arrangement of the faults.
To cope with this issue, we propose to apply a constant false-alarm rate (CFAR) method, borrowed from radar literature [23], which sets an adaptive threshold in order to keep the false-alarm probability constant. In particular, we adopt the cell-averaging CFAR (CA)-CFAR. By this method, the detection threshold for a specific array element under test is  determined by averaging over some neighbouring cells. More in detail, consider the vector c arranged in matrix form C ∈ C N x ×N y , with N x and N y being the number of elemental radiators along x and y, respectively, and N = N x N y . Then, the element at (t, s) corresponds to the n array element with n = (s − 1)N x + t, t ∈ (1, 2, · · · , N x ) and s ∈ (1, 2, · · · , N y ). VOLUME 10, 2022 The detection threshold for the pixel (t, s) of the reconstruction is evaluated by first averaging the reconstruction over neighbour cells belonging to a rectangular window of K = K 1 × K 2 elements, that is with I T being the rectangle window centered in (t, s) and sides K 1 + N guard1 and K 2 + N guard2 , and I guard the rectangle window centered in (t, s) and sides N guard1 and N guard2 . Y (t, s) actually represents an estimate of the background noise (interference). Eventually, the detection threshold is obtained by multiplying the estimated average by a scaling factor α T . The latter allows to control the false alarm probability by Hence, by setting PFA d , α T is evaluated from (19) and the detection threshold is finally set equal to A th = α T Y .
Note that in (18), the element under test and its immediate neighbours belonging to I guard are excluded by the averaging procedure. This is done in order to avoid biasing the estimation of the threshold towards excessively high values. In this regard, N guard1 and N guard2 must be properly chosen according to the psf function which describes how the fault reconstruct actually spreads over more than one pixel point. In [20], an analytical estimation in terms of the configuration's geometrical parameters is derived. This result can then be used to set the guard pixels as those for which the main-lobe of the psf is excluded by the averaging procedure. More in detail, it is shown that only when X 0 → ∞ and Y 0 → ∞, the main-lobe half-width of the psf tends to λ/2, otherwise it is larger and non-uniform across the array domain. In particular, for the cases considered in this manuscript, the psf estimate reported in [20] suggests to set N guard1 = N guard2 = 1.
The detection threshold resulting from (18) and (19) does not properly take into account the contributions of the reconstructions of the faults that can enter the averaging window. The result is that noise estimation is distorted and generally increases the threshold. This leads to a deviation between the set PFA d and the actual one. Accordingly, also K 1 and K 2 must be judiciously chosen. Indeed, K 1 and K 2 must be chosen so as to make the noise power estimation reliable across all of the array element positions. In principle, larger K 1 and K 2 imply better noise power estimation. However, this is not necessarily true when more than one fault is present. In fact, excessively large K 1 and K 2 can decrease the detection threshold reliability because other faults are included in the averaging due to the strong interference of the main lobes. After some trials, we set K 1 = K 2 = 10. Now we are ready to assess the detection performance. To this end, as a preliminary example, we just address the case reported in Fig. 3. In particular, since BTM and MM return similar results, we show only the detection concerning the MM method. The outcome of the detection is reported in Fig. 6, where SNR = 20dB and the detection threshold A th is set by fixing PFA d = 0.001. As can be seen, in both the cases, the faults are clearly and correctly detected and localized. Remarkably, this holds true also for the non-uniform sampling strategy. The case of Taylor distribution also shows many false positives (for both sampling methods).
These results show that the suggested data reduction strategy actually works. However, it is not conclusive since it refers to just one example. In order to estimate the probability of detection and of false alarm we run a Monte Carlo analysis by considering N trials = 1000 different realizations of noise and fault layouts. In particular, PFA and PD are evaluated as (20) and PD = 1 N trials N trials n=1 FD(n) N fault (21) where FP, FD and N fault denote the number of false positives, the number of detected faults and the number of faults, VOLUME 10, 2022  Fig. 4.

TABLE 2.
Detection results for the case of Taylor tapered excitation coefficients with SLL = −20dB andn = 4. The geometrical parameters are the same as in Fig. 5.
respectively. Different values of SNR and N faults are also considered for each sampling strategy, uniform and non uniform (nu). In particular, the results concerning the uniform distribution are shown in the Table 1, whereas Table 2 refers to the Taylor distribution. In both cases, the threshold is chosen by fixing PFA d equal to 0.001. As can be seen, uniform sampling generally returns results which are only very slightly better than the non-uniform sampling. Actually, they are really very close, and hence demonstrate the effectiveness of the non-uniform sampling strategy. Also, as expected, PD degrades when SNR reduces and/or the number of faults increases. For the Taylor distribution, as natural, this is a bit more evident. Nonetheless, in all the cases, the results keep to satisfying values, also considering that if antenna diagnostics is performed in controlled environments like an anechoic chamber, the SNR can be sufficiently high. Finally, it is noted that PFA decreases when the number of faults increases. This, at first glance, could appear contour-intuitive but it is actually a consequence of the chosen detection strategy. Indeed, when many faults are present (18) tends to return higher detection thresholds and hence lower PFA. As anticipated above, this is how fault reconstructions distort threshold setting.

VI. CONCLUSION
In this paper, the problem of planar array antenna diagnostics from near-field measurements has been addressed. In particular, the standard Back transformation and the Matrix methods have been considered in view of the new field sampling strategy proposed in [22]. The latter provides a non-uniform sampling grid that allows to reduce the number of measurement points compared to the standard λ/2 strategy, resulting in a decrease of the measurement acquisition time. In particular, for MM this positive impact becomes even more relevant because the reduction in the number of points also decreases the computational burden.
Since from the reconstructions it is quite difficult to identify the fault locations, a fault detection step has been added to the array diagnostic problem. Such a step introduces an adaptive threshold related to a reference current and relies to a cell-averaging CFAR (CA)-CFAR technique [23]. It allows to choose a detection threshold which keeps the false-alarm probability constant despite changes in the background noise.
Finally, numerical analysis has shown that the non-uniform sampling strategy yields results, in terms of probability of detection and of false alarm, which are comparable to the ones returned by the more standard λ/2 sampling, with much fewer data. Indeed, uniform sampling returns better results only for low SNR. This could be expected and it is the trade-off for reducing the number of data. This drawback can be partially remedied by increasing the number of data through a higher oversampling factor ν.