Accurate Reconstruction of the Radiation of Sparse Sources from a Small Set of Near-Field Measurements by Means of a Smooth-Weighted Norm for Cluster-Sparsity Problems

The aim of this contribution is to present an approach that allows to improve the quality of the reconstruction of the far-field from a small number of measured samples by means of sparse recovery using a relatively coarse grid for source positions (with sample spacing of the order of λ/8) compared to the grid usually required. In particular, the iterative method proposed employs a smooth-weighted constrained minimization, that guarantees a better probability of correct estimate of the sparse sources and an improved quality in the reconstruction, with a similar computational effort respect to the standard ℓ1 re-weighted minimization approach.


Introduction
The evaluation of the far-field pattern of radiating sources is one of the most important procedures for the correct design of modern communication systems. Moreover, the correct operation of antenna systems may require testing, validation, and diagnostics of the radiating equipment.
Near-field measurement systems allow to perform measurements in a controlled environment with great advantages in terms of accuracy and reliability of the test site [1]. However, data acquisition is generally time consuming.
For example, in planar NF systems a λ/2 spatial sampling step [2] is currently universally adopted. This gives a huge number of samples in case of large antennas. Consequently, a large effort has been devoted to reducing the amount of data that must be acquired. Important steps in this direction have been obtained using the theory of Number of Degrees of Freedom (NDF) of the electromagnetic field [3]. The characterization of the far-field of sources by mean of near-field measurements with a good accuracy requires, in general, a number of measurements that is slightly larger than the NDF of the electromagnetic field of a source of a certain size and dimension [3][4][5][6].
Recently, further approaches for field interpolation from a small set of measurements have been proposed. In particular, many of these approaches have been inspired, by the compressive sensing paradigm [7], that exploits the knowledge of the sparsity of the chosen representation of the source to reduce the number of needed measurements. Some of these approaches have been aimed to the estimation of the far-field pattern from a limited number of points [8][9][10][11][12][13][14][15][16][17]; other methods have been developed for the identification of the failures of the radiating elements of arrays [18][19][20][21][22]. A review of the NF measurement techniques is reported in [23]. The paper includes a description of the methods based on standard 1 minimization technique and a comparison among different approaches in terms of number of measurements required for an effective estimation of the field.
One of the features that is common to many of the aforementioned compressive sensing approaches is the representation of the source by using a grid of samples defined in the region in which we expect the sources/failures to be positioned. Although these approaches do perfectly fit the problem of failure identification (since the exact position of the radiating elements is information that can be achieved with an excellent accuracy), in other cases, as for example antennas covered by radome, the position of the source may not perfectly match the employed grid.
To overcome this problem a possible solution is to use a dense grids. This is equivalent to approximate the field radiated by Antenna Under Test by the field radiated by an array, that in the following will be called "Grid Array" (GA), whose radiating elements are in the intersection points of the grid. The correct operation of this kind of algorithms requires the use of very dense grids, with spacing as low as λ/100 [24] to achieve the wanted accuracy in field reconstruction or source identification, with a significant impact on the computational burden of the reconstruction procedure, that may become infeasible in case of large planar or volumetric arrays.
This suggests the use of a less dense grids. However, in this case the positions of the AUT sources generally do not coincide with elements of the GA, and the sources of the AUT may be approximated on the GA by "clusters" of sources around the AUT sources ( Figure 1). This is a critical problem in sparse recovery, since the number of measurements required by the sparse-retrieval algorithms are linearly related to the degree of the sparsity of the problem [7], i.e., to the number of non-zero values of the unknowns. Consequently, the presence of clusters of, for example, four non-zero elements in the GA associated to the single source in the AUT, increases the number of needed measurements of 4 times.
x y Figure 1. The unknown positions of radiating elements of the Antenna Under Test (AUT) are approximated by a grid of points, obtaining a "Grid Array" (GA) whose radiating elements are in the grid positions; the red points represent the position of a source of the AUT; if the position does not coincide with the position of one element of the GA, the AUT element excitation may be approximated in the GA by a cluster of radiating elements (red squares); the effect is an increase in the number of non-null radiating elements of the GA, i.e., of the order of sparsity of the retrieval problem; this causes an increase in the number of measurements required to solve the sparse problem, that is proportional to the number of non-null elements of the problem.
The aim of this contribution is to propose an algorithm that allows a good quality of the reconstruction of the far-field using a relatively coarse grid (with sample spacing of the order of λ/8). The approach followed in this paper is to introduce a proper weighted norm that allows to identify the number of clusters in the GA grid. This allows to consider the presence of the clusters instead of simply the non-zero elements in the vector of the unknowns, improving the accuracy of the reconstruction.
The paper is organized in the following way: in Section 2 we provide the problem formulation, we discuss one of the existing approaches to the problem and we propose a novel method. In Section 3, we analyze in detail a specific problem. In Section 4, we provide a number of statistical tests to validate the proposed procedure. The numerical examples discussed in Sections 3 and 4 regard sparse sources; it is understood that these sources can be "true" radiating elements, or "equivalent" radiating elements in case of application of differential techniques for array diagnosis [9].

Problem Formulation
In this paper we will consider harmonics sources working at frequency f . The corresponding free-space wavelength and wave-number will be λ and β, respectively.
Let us consider as Antenna Under Test (AUT) a sparse array consisting of S radiating elements with unknown positions and excitations. For the sake of simplicity, in the following section we will consider point-like radiating sources. The S elements are placed in unknown positions (x s , y s , 0) (s ∈ [1, S]), in a rectangular area Ξ in the range x ∈ [x min , x max ] and y ∈ [y min , y max ], of the (x, y, 0) plane as in Figure 2. The field radiated by the AUT sources is measured in M points having coordinates (x m , y m , z m ) in radiative near-field area of the antenna [2]; the measurements are supposed to be affected by additive white Gaussian noise.
Accordingly, the field in the measurement positions is given by: whereinb ∈ C M is a vector collecting the value of the field measured in the M positions, s ∈ C S is a vector collecting the excitations of the S sources, n ∈ C M is a complex AWGN vector modelling the presence of noise, andÂ ∈ C M×S is the radiation matrix that relates the fields in the observation points to the excitations, whose elements are: and R m,s is the euclidean distance between the s-th source and the m-th observation point.
The goal is the reconstruction of the overall far-field radiated by the AUT from the knowledge of the field in the M measurement points. Different choices are obviously possible with inessential changes to the overall discussion.
A possible approach to solve this task consists in retrieving the excitation coefficients of the AUT from the M measured data, and to estimate the far-field (or the field in any point of the radiative area of the array) from the retrieved excitations. This method is effective and simple [25] but requires the exact knowledge of the positions of the AUT radiating elements.
As briefly discussed in the Introduction, in order to overcome this problem it is possible is to introduce an equivalent "Grid Array" (GA) considering a dense regular grid of N = N x × N y radiating sources with positions (x n ,ỹ n , 0) and inter-source spacing δ << λ/2 that sample Ξ with a sufficient approximation. The excitation coefficients of the grid array will be collected in a complex matrix The field radiated by the GA and calculated in the same positions where the radiated field of the AUT is measured are collected in a complex vector b ∈ C M .
According to the above model, the field radiated by the GA in a measurement point is wherein "vec" is a linear operator that performs the stacking of the columns of a matrix into a single column (in this case transforming a C N x ×N y into a C N column vector), A ∈ C M×N is the radiation matrix that relates the fields in the points b ∈ C M to the excitations vec(X) ∈ C N , whose elements are being R m,n is the euclidean distance between the n−th source point on the rectangular grid and the m-th observation point.
The main idea is to estimate an equivalent GA source that emulates at best the field of the unknown AUT: where is a threshold whose value depends on the noise n m affecting the measurementsb. Estimation of X from the M measurements is straightforward if M ≥ N, and could be achieved by means of pseudo-inversion [26]. However, in our case the number N of elements of the GA is much larger than the number M of the measurements points. Moreover, even if we considered a sufficiently large number of measurements, the inversion of the linear operator A would be affected by a severe ill-conditioning because of the limited NDF [3].
Actually, there is an infinite number of different GA sources that are able to emulate the same measurement vectorb and we are not guaranteed that when finding one of these GA their far-field would be the same of the real AUT. Proper additional conditions to regularize the inversion are indeed needed.

The Standard Re-Weighted Approach
As preliminary step, we recall the re-weighted sparse minimization algorithm [27]. Let us consider a dense regular grid of N = N x × N y sources. The reconstruction of the coefficients of the GA can be achieved solving the following problem, involving the so called "zero norm", i.e., the number of non-null elements of a vector or matrix.
As is well known, this minimization, because of the non-convexity of the 0 norm, is an NP problem. One of the currently adopted solutions to solve this minimization is to recast the problem as a re-weighted 1 norm minimization [27][28][29], that may be obtained with the following iterative constrained minimization: wherein X • W 1 is the Hadamard entry-wise product of the two matrices, * p with p ≥ 1 is the p norm of a matrix/vector, equal to the p−th root of the sum of the amplitudes of its entries to the p−th power, W is a matrix of weights whose entries X i,j are usually defined as where η a small constant, of the order of the minimum excitation we want to reconstruct [27]; as weight for the first iterations, because of the lack of a previous solution, a constant matrix is employed. The minimization (8) and (9) can be repeated for a small number of times, usually less than ten times, until the convergence is reached (i.e., the solution found does not change significantly between successive iterations). Because of the sparsening and regularizing effects of the weighted 1 minimization the solution found is usually sparse (i.e., few entries of the found matrix X are non-null) and may match the position of the actual unknown S sources. It has to be underlined that the aforementioned procedure may fail when the ratio M/S is not sufficiently high and/or when the amplitude of the noise is high [27].
Once the source has been reconstructed, it is possible to calculate its pattern, f.i., its far-field pattern, thus achieving the wanted reconstruction.
As discussed in the Introduction, the previously described approach usually requires very dense GA, with spacing of the order of λ/100, that limits the dimension of the AUT that can be tested with this technique.
This issue suggests the use of less dense grids, f.i., choosing a λ/8 step. However, in this case, the positions of the S sparse sources may not coincide with elements of the GA, and a good approximation of the S sparse source could be found by means of "clusters" of non-zero excitations in the GA [30,31]. The standard weighting vector is not proficient in dealing with such excitation clusters.

The Smooth Re-Weighted Approach
Inspired by the excellent results obtained for the antenna array synthesis with the smooth re-weighted L1 minimization, introduced in [31,32], a similar approach will be now proposed for the solution of the source excitation. The main idea is to properly construct the weighting matrix so that we can deal correctly with the clusters that could approximate the real pattern with a coarse grid; in the following, we will call "cluster-sparse" a grid excitation that is non-null apart from some some 2 × 2 cluster sub-matrices, and each cluster has a distance of at least 1 row and/or 1 column from the others.
The proposed novel approach works in this way: in the first iteration the same convex problem of (8) and (9) will be solved, but in the successive iterations a different task would be faced, in which the smooth-weighted norm is used as a constraint and not as the term to be minimized: subject to where W D is the "smooth" weight matrix, and ν D is a proper threshold, whose meaning will be clarified in the following; the smooth weight matrix is calculated by first defining the smooth excitation matrix X D = abs(X) * D (13) where "abs" takes the absolute values of the entries of matrix X, * performs a 2D matrix convolution that outputs only the internal part of the convolution (of the same size of X) with then the entriesW D i,j of the smooth weight vector are calculated as: where τ is a small constant that prevents (15) to become too large when small or null values of X D i,j are in order. The rationale of the use of the weight in (15), that does not depend on the value of the excitation "point-by-point", but is "smooth" and non-local, is due to the fact that when the re-weighted minimization has converged, if the excitation matrix X is cluster-sparse, then where N D is the number of clusters in the matrix X. A demonstration for the linear case, whose 2D case extension is straightforward, is provided in [31]; a simple numerical example of this feature is discussed in Appendix A. Obviously, clusters larger than 2 × 2 can be managed by this approach by considering a larger D matrix in (14).
Let us now discuss the convex iteration (11) and (12). Differently from the approach described by (8) and (9), that was used in previous works as [23,32], the quantity to be minimized is now the difference of the measured radiated field with respect to the field radiated by the dense source, and the weighted 1 norm is now used as a constraint that, as previously discussed, limits the number of clusters in the excitations. At each iteration, before solving (11) and (12) we can calculate the value of ν D = X • W D 1 using the previously calculated excitation.
After a number of iterations is performed, and an estimate of the source has been reconstructed, it is possible to calculate also in this case the radiated field, thus achieving the wanted sparse source pattern retrieval.

Numerical Examples
In this section, we will analyze an example, in order to better understand the proposed method.
The region Ξ will be a square in the (x, y) plane, with about 8λ side; The regular grid of N = 4096 points on this plane will have a spacing of δ = λ/8. As sparse sources we will consider 5 isotropic sources, with their coordinates in the vertices of a regular pentagon, with a radius of 3.5λ. The pentagon is rotated of 18 • in the (x, y) plane, and its excitations have been chosen equal in amplitude, but with a random phase; the coordinates and phases of the sources are provided in Table 1. The data are collected on M = 25 measurement points on a spherical surface having a radius equal to 8λ using a spiral scanning [33,34]. The coordinates of the measurement points are provided in Table 2. A sketch of the system is provided in Figure 3. It is worth noting that in case of non-sparse sources a number of 4π3.5 2 150 measurement data points are required [4].
We will now apply the two sparse source identification algorithms discussed in the previous section (the "standard" one and the "smooth" one). Both the algorithms will be executed for a fixed number of iterations N i = 10 on the same aforementioned dense grid; both algorithms would process the same measurement vector (1), that has been simulated with an SNR = 40 dB, calculated as the ratio of the norm of the received field and the norm of the noise. The threshold will be 2 times the value of the norm of the noise, and τ = η = max(X i,j )/100. For the computation of the minimization problems we have used the CVX package for Matlab [35,36]   In Figure 4, we provide the imagemaps of the matrices representing the reconstructed sources in the two cases. As we can see, both algorithms are able to identify the approximate position of the 5 sources, but the reconstruction with the proposed algorithm seems to be more accurate, since the sources are reconstructed as "clusters" instead of single points on the grid, thanks to the "smoothing". It has also to be noticed that there are also some "sparse" artifacts in the reconstruction with the standard algorithm, that are not present in the proposed approach. Let us now look at the reconstructed fields: in Figure 5 we can compare the "true" far-field with the far-field calculated from the two reconstructed sparse sources. The reconstruction obtained in the case of the standard algorithm is of sufficient quality, but the reconstruction with the smooth algorithm is by far superior, and almost indistinguishable from the "true" field. To quantify the difference of the two fields, in Figure 6 we provide an imagemap of the difference of the reconstructed fields and the true field. The plot confirms the superiority of the "smooth" approach with respect to the standard one, but to achieve a precise quantification of the reconstruction error we calculate the relative error as: In the considered case, we get the values of ∆ FF,standard = 0.306 and ∆ FF,smooth = 0.033.  It is worth mentioning that the computational complexity of both algorithm is equivalent, since they are executed with the same GA (the calculation of the smooth weight is negligible from a computational point of view). It must be also noted that different choices could have been completed for the values of , τ, and η, but according to our tests, the final results would not change significantly; similarly, we have tested other kinds of field sampling rules for the positioning of the measurement points, apart from the employed spiral ones, finding similar results. Obviously, the presented case is only a single test, from which we could not claim the superiority of the proposed approach. For this reason, in the next section we will provide a detailed statistical analysis.

Statistical Analysis
In this section, we will compare the proposed "smooth" weighted algorithm for the sparse source field reconstruction with respect to the standard approach.

Statistical Analysis with 5 Random Sources
In this preliminary test, we will consider the same setup of the example in the previous section, but instead of considering a single set of random sparse sources, we will analyze P = 1000 random sets of 5 isotropic sources placed in the Ξ region; the position of the isotropic sources will be obtained with a uniform distribution of probability, but sets with inter-source distances lower than λ/2 will be discarded.
In Figure 7, we can see the comparison of cumulative density function of the reconstruction error (17) for the standard algorithm and the proposed algorithm, ran with the same noise realization (with a SNR of 40 dB): the "smooth" algorithm shows a much smaller reconstruction error on the far-field.
It is also worth observing that the CDF for the smooth algorithm shows a "jump" in the error rate for a value of the CDF of 0.96. This behavior is not present in the distribution of error of the standard algorithm that is much more "regular". When these jumps occur, the reconstruction is completely messed and allow us to analyze the correct detection rate: the proposed algorithm allows the detection of up to 5 sources in 96% of the cases.

Analysis with Variable Number of Sources
In this second test, we will analyze what happens with a variable number of sources S, varying from 1 to 10, always considering the same SNR of 40 dB. As we can see the difference in the reconstruction error of the two approaches confirms the validity of the proposed approach.
The comparison, provided in Figure 8 confirms the result that was obtained in the previous example; in particular while the standard algorithm shows a smooth behavior of the CDFs, the proposed approach shows a "jump" from S = 5 to S = 9; the detection rate is greater than 99% up to four sparse sources, and is less than 1% for ten or more sources. Table 3 collects the feasible reconstruction rates (FRR), i.e., the values of the CDF for which the jump occurs.
These results confirm that the proposed approach is able to detect a higher number of sparse sources with the same number of points, since the weighting approach employed deals with clusters such as a single entity, differently from the standard approach for which the presence of clusters may increase the level of sparsity of the problem from S to αS, where α > 1 is approximately equal to the average dimension of the clusters used for the approximation of the AUT in the GA. The standard reconstruction algorithm would require a higher number of measurement points M or a denser GA to achieve an accuracy comparable to the proposed approach.

Analysis with a Variable Noise
In this last example, we will compare the two algorithm considering a set of S = 5 sources and a variable noise. In particular, we will focus on a variation of the SNR from 30 dB to 50 dB.
In Figure 9, we show the behavior of the two approaches; as expected, the good performance of the proposed algorithm with respect to the standard one are confirmed. It is also possible to see that the variation of the SNR marginally changes the detection rate, and affects only the values of the reconstruction error ∆ FF . Figure 9. Cumulative density function of the reconstruction error for the random sets of S = 5 sparse sources in the Ξ region for a variable SNR. Solid lines: "standard" approach; dashed lines: "smooth" approach.

Conclusions
In this contribution, we have analyzed the problem of the reconstruction of the field radiated by sparse sources, whose radiating elements positions are not known, using a reduced number of measured samples.
The algorithm proposed in this paper is based on the estimation of the sparse source by using the smooth weighted 1 norm. The approach proposed allows to effectively exploit relatively coarse grids used for the approximation of the sparse sources, allowing better reconstruction performance with respect to the standard re-weighted 1 minimization algorithm.
Numerical simulations confirm that using the same number of measurement points and the same computational effort, it is possible to reconstruct the far-field of the AUT with a significantly better accuracy with respect to the standard approach.
The analysis of this and the previous section has been performed for isotropic sources. This choice has been completed for limiting the computational effort required for the statistical analysis, and for simplifying the comprehension and presentation of the results, that can be easily repeated by the interested reader.
As a future development, we are planning to apply the proposed approach also to the detection of failures of linear, planar and conformal arrays, and also to use similar techniques for the rapid testing of arrays. Finally, we plan to perform further comparisons with other approaches, as well as considering more realistic models to test the proposed approach. Funding: This paper has been supported by the MIUR program "Dipartimenti di Eccellenza 2018-2022".

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Weighted Norm for Cluster-Sparsity Problems
The aim of this Appendix is to clarify the properties of the weighted norm proposed in this paper. For sake of clarity, the characteristics of the norm will be discussed using a numerical example.
In this section, we consider a simple excitation matrix X. As preliminary step, before we discuss the proposed weighted norm, we will calculate the weight vector with the standard approach.
Let us now consider As we can see, in the standard "penalty" introduced by the weight vector is inversely proportional to the values in X: if we calculate X • W 1 we would get 9.9363 ≈ 10, thus verifying that the re-weighted norm asymptotically tends to approximate the 0 norm of the X, i.e., the number of non-null entries of X. The fact that the excitation are organized into clusters is absolutely not taken into account: the 1 weighted norm returns a measure of the sparsity that is much larger with respect to the number of clusters.
If we followed the smooth weighted approach with η = 0.01 we would get: The weight corresponding to each one of the elements of any cluster is inversely proportional to the sum of the excitations of that cluster. If we calculate the 1 smooth weighted norm, namely X • W D 1 , we would get 4: the measure of the sparsity induced by the smooth weight is the number of clusters of dimension non greater than 2 × 2 in X-as we wanted. Using the smooth weighted norm in the minimization, or as constraint, allow us to correctly represent the sparsity of the source when the sparse sources are approximated as clusters of sources.