On the Sampling of the Fresnel Field Intensity over a Full Angular Sector

: In this article, the question of how to sample the square amplitude of the radiated ﬁeld in the framework of phaseless antenna diagnostics is addressed. In particular, the goal of the article is to ﬁnd a discretization scheme that exploits a non-redundant number of samples and returns a discrete model whose mathematical properties are similar to those of the continuous one. To this end, at ﬁrst, the lifting technique is used to obtain a linear representation of the square amplitude of the radiated ﬁeld. Later, a discretization scheme based on the Shannon sampling theorem is exploited to discretize the continuous model. More in detail, the kernel of the related eigenvalue problem is ﬁrst recast as the Fourier transform of a window function, and after, it is evaluated. Finally, the sampling theory approach is applied to obtain a discrete model whose singular values approximate all the relevant singular values of the continuous linear model. The study refers to a strip source whose square magnitude of the radiated ﬁeld is observed in the Fresnel zone over a 2D observation domain.


Introduction
In the framework of inverse problems in electromagnetics [1][2][3][4][5][6][7], the inverse source problem plays an important role. The latter is a classical problem of the electromagnetic literature [8][9][10][11][12] which consists in recovering the source current from the observation of its radiated field. From the mathematical point of view, the inverse source problem requires to invert a linear integral operator, called a radiation operator, that relates the density current J with the radiated electric field E.
Despite the performances of electronic instrumentation are continuously improved, in some contexts the measure of phase information may be not accurate. Hence, it is worth addressing the problem of reconstructing the source current from amplitude-only data of the radiated field [13][14][15][16][17][18][19][20]. Such a problem falls into the realm of phase retrieval, and it can be formulated as recovering the density current J from the quadratic model: where T stands for the radiation operator. One of the most popular algorithms to find a solution consists in the minimization of a cost functional obtained by the least-squares problem associated to Equation (1) (see [21,22]). Since the functional to be minimized is quartic, such a method may suffer from the local minima problem.
Different studies have shown that a crucial parameter in the analysis and cure of local minima is played by the dimension of data space M, which represents the number of independent functions that allow expressing the data with a given degree of accuracy [23]. According to such studies [21,24,25], once the number of unknowns N has been fixed, the functional to be minimized becomes more and more similar to a convex one when the dimension of data space M increases. Hence, it will exist a particular value of the dimension of data space such that the functional is free from traps. It is worth noting that the question of quantifying the minimum value of M that ensures the convergence of the least-squares approach is still an open point, since such a value has been found only for some particular mathematical models [24][25][26][27].
In light of the previous discussion, the first aim of this paper is that of evaluating the dimension of data space M and establishing how such quantity is related to the geometrical parameters of the configuration.
Once the dimension of data space will be evaluated, the second goal of this article is that of providing an efficient discretization of the continuous model shown in Equation (1). This study is motivated by the fact that the square amplitude of the radiated field cannot be continuously collected over the observation domain; accordingly, it is of great practical interest to find a discretization scheme that: • requires a non-redundant number of measurements; • allows to obtain a discrete model that shares the same mathematical properties of the continuous one.
Such a discretization scheme would bring a series of advantages. On one hand, it would reduce the acquisition time to measure the data, and the memory requirements to store them on the other, it would decrease the computational burden to perform the inversion.
Regarding the sampling of the radiated field E, standard sampling schemes for the case of planar, cylindrical, and spherical scanning are described respectively in [28][29][30]. Despite this, such strategies are not efficient since they require to acquire a number of samples that may be significantly higher than the dimension of data space. Indeed, especially at microwave and millimeter frequencies, the number of required measurements may become extremely high.
In order to reduce the number of measurements, during the years a different sampling scheme (mostly based on a spatially varying spacing) have been developed. Some of them exploit the concept of local bandwidth [31][32][33], others use an adaptive procedure [34][35][36], still others optimize some metrics related to the singular values decomposition of the radiation operator [37][38][39][40][41][42]. In addition to those just mentioned, further sampling techniques used in antenna applications are described in [43][44][45][46][47].
A particularly interesting discretization scheme is the sampling theory approach developed in [45,46]. The latter, under the hypothesis that the operator TT † is convolution and band-limited (with T † denoting the adjoint operator of T ), provides not only an efficient sampling of the radiated field (or in other words the grid of points where the data must be collected) but also a strategy to discretize the linear model T J = E.
However, the sampling approach is suitable only for linear problems; hence, it cannot be directly applied in the case of phaseless measurements. With the aim to overcome this issue, the lifting technique [23] can be exploited to recast the quadratic model in (1) as a linear one. Once a linear representation of |E| 2 has been obtained, the sampling theory approach can be fruitfully employed to efficiently discretize the continuous model.
Let us stress that the direct applicability of the sampling theory approach is limited to the cases where the operator involved in the correspondent eigenvalue problem is convolution and its kernel is a band-limited function. Hence, when such an operator does not fulfill these properties, the sampling theory approach cannot be directly applied. In such cases, indeed, it can be employed only after a warping transformation has been exploited to recast the kernel of the operator involved in the correspondent eigenvalue problem as a band-limited function of a difference type [48][49][50].
The study developed in this paper will be done with reference to a 2D geometry consisting of a strip electric current observed on an extended observation domain in Fresnel zone.
The paper is organized in the following manner. In Section 2, the geometry of the problem and the correspondent radiation operator are introduced. In Section 3, at first, the lifting operator is defined; after, its analytical properties are studied. In Section 4, the dimension of data space is estimated. In Section 5, a discrete model that approximates the relevant singular values of the lifting operator is provided. A section of conclusions follows.

Geometry of the Problem
Let us consider a 2D scalar geometry consisting of a source whose density current J(x) = J(x)î y is supported on the set [−a, a] of the x-axis and directed along the y−axis (see Figure 1). Such current radiates an electric field E(r, θ) = E(r, θ)î y in a homogeneous medium whose wavenumber β is given by β = 2π λ with λ denoting the wavelength. The square magnitude of the electric field |E| 2 is observed in the Fresnel region on a 2D domain that extends along the polar coordinates (r, θ) on the set [r min , For the considered geometry, the radiated electric field E can by expressed by the equation where u = sin(θ), and the radiation operator T is such that with L 2 (SD) denoting the space of square-integrable functions on SD = [−a, a], and L 4 (OD) indicates the space of functions whose amplitude to the fourth power is integrable Under the paraxial Fresnel approximation, Equation (2) can be explicitly written as

Study of the Lifting Operator
The first aim of this section is that of providing a linear representation of the square magnitude of the radiated field |E(r, u)| 2 . The introduction of a linear model allows addressing two important tasks. In particular,

1.
it enables to estimate the dimension of data space by exploiting a singular values approach; 2. it allows finding a discrete model that shares the same singular values of the continuous one.
With the aim of obtaining a linear representation of |E(r, u)| 2 , the lifting technique can be exploited. The latter adopts a redefinition of the unknown space to recast the quadratic model |T J| 2 = |E| 2 as a linear one. In particular, a linear model is obtained by rewriting Equation (1) in the form below and, by considering as unknown the function In such a way, the integral operator A which links the unknown function F(x, x) with the data function |E(r, u)| 2 is linear. The latter is known in the literature as the lifting operator and it is defined as where Accordingly, the square magnitude of the radiated electric field can be expressed as Once that linear representation of |E| 2 has been achieved, our goal is that of addressing the tasks mentioned in points 1 and 2. To this end, it is necessary to introduce the adjoint operator A † and to explicit the auxiliary operator AA † . Indeed, as well known, the singular values of A are equal to the square root of the eigenvalues of AA † [51].
The adjoint operator A † is usually defined as where (·) denotes the function of the variables (r, u) on which the adjoint operator acts. In Appendix A, the integral operator AA † is introduced and a closed-form expression of its kernel is provided.
Differently from what is done in Appendix A, here, a weighted adjoint operator A w is exploited. The latter can be expressed as where w(x, x) represents a weight function. On the one hand, the use of a weighted adjoint will allow simplifying the mathematical discussion and on the other, it will bring a change of the singular values of the lifting operator. Despite this, such changes will affect only the dynamics of the singular values but not the critical index at which they become negligible. Hence, for our purposes, the effect brought by the weight function can be neglected.
Anyway, the main difference with what is done in Appendix A is that, here, the kernel of the operator AA † w will be expressed as a double integral and not as the square amplitude of a one-dimensional integral. This new representation of the kernel of AA † w will suggest to introduce a change of variables or better a warping transformation that allows rewriting it like the Fourier transform of a constant function defined on a two-dimensional finite domain. On the basis of such discussion, it follows that: • the eigenvalues of AA † w can be computed in closed-form by resorting to the Slepian Pollak theory; • the kernel of AA † w will be not only convolution but also band-limited, hence, the sampling theory approach will be fruitful to find a discretization of the correspondent eigenvalue problem.
Note that the idea of exploiting a warping transformation with the aim to recast the kernel of an operator like the Fourier transform of a window function has been already exploited in other recent works [48][49][50].
In order to formalize the above discussion from a mathematical point of view, let us explicit the operator AA † w . The latter is given by where: Now, let us divide the integration domain S = [−a, a] × [−a, a] as made below Since S 2 is a null set with respect to the Lebesgue measure, the kernel H(r, r o , u, u o ) can be computed by performing the integration only on S 1 . This means that With the aim to express the kernel through a Fourier Transform relation, it is possible to introduce the following transformation The latter is injective and continuously differentiable on S 1 and, it allows rewriting the kernel of AA † w in the following form where • Σ 1 denotes the new domain in which the original integration domain S 1 is mapped by trasformation (15), is the Jacobian determinant of the transformation.
Let us remark that, despite ∂(X 1 ,X 2 ) being singular for X 1 = 0, such point does not belong to the integration domain Σ 1 . Hence, the integrand in (16) is free from singularities.
At this point, let us choose the weight function w(X 1 , X 2 ). By fixing it results that Equation (18) shows that the kernel of AA † w can be seen as the Fourier transform of a constant function defined on the set Σ 1 . The shape of the set Σ 1 is depicted in Figure 2 . As seen from Figure 2, the integration domain Σ 1 is not rectangular. However, since the area of Σ 1 does not differ significantly from the area of the smallest rectangular that encloses Σ 1 , H(r, r o , u, u o ) can be evaluated by integrating on the smallest rectangle that includes the set Σ 1 . The latter is made up by all the points (X 1 , X 2 ) belonging to the rectangular set [−2a, 2a] × [−a 2 /r max , a 2 /r max ].
By performing the integration on such domain, it results that Accordingly, the operator AA † w can be expressed as From (20), it is evident that the operator AA † w becomes more similar to a convolution operator if s = r max /r is set. In fact, by doing this, it follows that (21) where s min = 1, s max = r max /r min . The last equation provides an approximated form of the operator AA † w which will be used in next sections to estimate the dimension of data space, and to find a discrete model that shares the same mathematical properties of the continuous one.

Estimation of the Dimension of Data
In this section, an analytical evaluation of the dimension of data space is provided. Different methods can be exploited to achieve this task, here, an SVD-based approach that allows estimating the dimension of data space by counting the number of significant singular values of the lifting operator A is exploited [52].
As is well known, the singular values of A, {σ m }, are related to the eigenvalues of AA † , {λ m }, by the equation However, only the eigenvalues of AA † w are known in closed-form. For such a reason, at first, it will be solved by the eigenvalue problem where λ m and v m are respectively the eigenvalues, and the eigenfunctions of AA † w . After, through a numerical analysis, it will be checked that the number of relevant eigenvalues of AA † w is essentially the same as the number of relevant eigenvalues of AA † .
To evaluate the eigenvalues of AA † w , let us explicit Equation (23) as below 8a 3 By fixingṽ Equation (24) can be recast as The eigenvalues of (26) were computed in closed-form in [53], and they are given by where λ it results that the eigenvalues of the problem (26) and, consequently, the eigenvalues of the operator AA † w are significant until the index Let us recall that the kernel of the operator AA † w has been computed by integrating on the smallest rectangular that encloses Σ 1 . For such a reason, the scalar M is not exactly equal to the number of relevant eigenvalues of AA † w but it represents an upper bound. In the first part of this section, the eigenvalues behavior of AA † w has been found, however, the actual singular values of A are equal to the square root of the eigenvalues of AA † . This means that only an approximation of the singular values of A given by the square root of the eigenvalues of AA † w is known. Now, by means of numerical experiments, it is checked that actual singular values of A and its approximated version computed starting from the eigenvalues of AA † w become negligible at the same index.
As a test case, a configuration in which a = 10λ, u max = 0.5, r min = 25λ (s max = 4), r max = 100λ and (s min = 1) is considered. With reference to such a configuration, in Figures 3 and 4 the actual singular values of A and their approximation in the sense clarified above are sketched in natural scale and in dB. In particular, the blue, red, and black diagrams sketch respectively: • the square root of the eigenvalues of the approximated version of AA † w provided by (21), • the square root of the eigenvalues of AA † w , • the square root of the eigenvalues of AA † .
As can be seen from Figures 3 and 4, the square root of the eigenvalues of Operator (21) exhibits a multi-step behavior and it is relevant until index M = M u M s = 164.
However, our aim is to predict the critical index at which the actual singular values of A become negligible. By observing the behavior of the actual singular values (black diagram in Figures 3 and 4), it is evident that the singular values beyond index M = 164 are surely negligible while those before are almost all significant if the noise level is not so high. This implies that the use of the weighted adjoint changes only the dynamics of the singular values but not the critical index at which they become negligible. For such a reason, it is possible to state that M = M u M s is an upper bound for the dimension of data space that is very close to its actual value.

Sampling Approach
In this section, the problem of how to discretize the continuous model AF = |E| 2 in such a way that the discrete model approximates the most important singular values of the continuous model is addressed.
Since the singular values of A are the square root of the eigenvalues of AA † , this task can be recast as finding a proper discretization of the eigenvalue problem where v m stands for the eigenfunctions of AA † . For the sake of simplicity, instead of discretizing Equation (30), a discretization of the eigenvalue problem AA † w v n = λ m v n that approximates all the relevant eigenvalues of AA † w will be found. Afterwards, by means of some numerical simulations, it will be checked that the discretization scheme adopted for the operator AA † w works also for AA † . As shown in Section 4, by fixingṽ m (s, u) = v m (s, u)/s, the eigenvalue problem AA † w v m = λ m v m can be recast as The continuous model in (31) involves a convolution operator with a band-limited kernel, accordingly, it can be discretized by exploiting the sampling theory approach developed in [45,46].
A key role in the discretization process is played by the sampling frequency. Since the bandwidth of the kernel function with respect to the variables u o and s o is given by W u = 2βa and W s = βa 2 /2r max , the sampling steps ∆u and ∆s to be used in the discretization process must satisfy the conditions where χ ≥ 1 is an eventually oversampling factor. The application of the sampling theory approach provides the following discrete model where • L is a matrix made up by the scalars L pqij given by •ṽ m is the vector whose elements are the samples of the eigenfunctionṽ m (u, s).
In Appendix B, all the mathematical steps that allow passing from the continuous Model (31) to the discrete Model (34) are shown.
According to [45,46], apart for the truncation error in the sampling expansion, the matrix L exhibits the same eigenvalues of the continuous model.
Let us point out that the dimension of matrix L and of the vectorṽ m are equal to the number of sampling functionsM used in the expansion of the kernel function and of the eigenfunctions. Since, in our case, only the sampling functions corresponding to the samples falling into the set [−u max , u max ] × [s min , s max ] has been considered, it results that M =M uMs (36) whereM As can be seen from the Equations (36) However, as specified before, our aim is that of finding an efficient discretization of AA † . For such reason, it is verified if the sampling steps ∆u and ∆s used in the discretization of AA † w work also for AA † . Figure 5 shows the eigenvalues of AA † , and those of the correspondent discrete models obtained by using as sampling steps (∆u = π/W u , ∆s = π/W s ) and (∆u = π/(1.2 W u ), ∆s = π/(1.2 W s )). As can be seen from Figure 5, a sampling frequency exactly equal to the Nyquist rate is already sufficient to approximate the most important eigenvalues of AA † ; instead, a sampling frequency higher than the Nyquist rate allows for also approximating the eigenvalues whose value is lower. It is interesting to highlight that the sampling approach returns also the locations where the phaseless data (the square amplitude of the electric field) must be collected. is mapped into a spatially varying sampling in the variables (r o , θ o ). Such behavior is confirmed also by Figure 6 which sketches the grid of the optimal sampling points in the case χ = 1. Figure 6. Grid of the optimal sampling points in the case χ = 1.

Conclusions
In this article, two issues which fall into the realm of phaseless inverse source problem were addressed. The first one was that of estimating analytically the dimension of data space by means of the singular value decomposition of the lifting operator. The second one was that of introducing a strategy to collect the phaseless data. In particular, in regards to the first point, a good upper bound to the dimension of data space is provided and the relation between the dimension of data space, the frequency, and the geometrical parameters of the configuration found. In regards to the second point, a discretization strategy that allows obtaining a discrete model whose singular values approximate very well the most important singular values of the lifting operator was introduced.
The analysis was performed in the case of a 2D scalar geometry consisting of a strip electric current observed on a two dimensional observation domain in Fresnel zone.  (A13) The linear system (A13) can be easily recast in the matrix form Lṽ m = λ mṽm .