Phaseless Characterization of Compact Antenna Test Range via Improved Alternating Projection Algorithm

: A practical compact antenna test range (CATR) requires good quiet zone quality for antenna characterization. This paper addresses the phase proﬁle of the CATR quiet zone from the known intensity pattern of spatial domain and Fourier domain based on a combined alternating projection algorithm. The proposed algorithm is composed of Gerchberg–Saxton (GS) and Hybrid Input–Output (HIO) algorithms and the two algorithms with spatial phase perturbation (SPP) work collaboratively or independently under predesigned conditions. It is observed that the algorithm with random initial phase guess can always converge to an optimal solution by performing a series of hierarchical optimizations of the problem. The numerical results are in good agreement with simulated results in different frequency bands, overcoming the phase retrieval limitation of local convergence in the iterative process. Furthermore, to validate the effectiveness and robustness of the proposed procedure, the related discussions corresponding to different sampling areas in Fourier domain and different signal to noise ratios (SNRs) are given.


Introduction
With CATRs applied in the millimeter (mm) or sub-millimeter (submm) wave bands widely [1][2][3][4][5][6][7], increasing importance has been attached to the quiet zone disturbance assessment. In order to ensure effective CATR application, ±0.5 dB/±5°maximum amplitude/phase deviations from a constant value are allowed in the design and the implementation of a CATR. As is well known, due to variation of the temperature, probe positioning inaccuracy, and cable flexing, phase data acquisition becomes difficult, especially at mmand submm-wave bands. Consequently, the phase retrieval (PR) of the field through amplitude-only has important research value. Additionally, the PR problem is of great concern in many other fields such as electron microscopy [8], astronomy [9], crystallography [10], lithography [11], optical imaging [12], and antenna characterization [13].
For simplicity, a very broad class of PR problems can be expressed as follows, T[x(t)] = |y(w)|e j(w) . t and w denote coordinate vectors of object and observation domains, respectively. T is an operator, and x(t) is an unknown signal. |y(w)| is the amplitude of observation domain. Then, we find x(t) from |y(w)| and some a priori information. The alternating projection method known as an iteratively plane-to-plane backpropagation algorithm has been widely used. This class of methods originated from the classical Gerchberg-Saxton (GS) [14] algorithm where T is a Fourier transform operator. For the GS algorithm, the sets of attributes defined in the two domains are nonconvex, resulting in easily subject subjecting to trap into false solutions. For this reason, Fienup introduced a broad framework for iterative algorithms based on the GS algorithm in a series of papers, in which three main classes of algorithms, error-reduction (ER), basic input-output (BIO) as well as hybrid input-output (HIO), were proposed [15][16][17]. In addition, many modified methods based on the above alternating projection algorithms have been proposed for solving antenna phaseless measurement problems [18][19][20][21][22][23][24]. These methods retrieve phase information through two sets of amplitude. However, they are very sensitive to the initial phase guess. Besides, the size of each plane needs to be much larger than that of the required recovery aperture and the sampling rate should be at least λ/4.
Over the past few years, PhaseLift in a new PR framework provided by Candès [25,26] has been recently developed and has aroused widespread interest. Different from the alternating projection method for solving non-convex problems directly, the methods under the new framework rely on the relaxation of the original nonconvex problem into a convex one to get the global optimal solution. It is shown in PhaseLift [25,26] that lifting the PR problem to a higher dimensional and using semidefinite relaxation can transform the PR problem into a convex quadratic programming problem. When the measured signals are sparse signals, the stability of PhaseLift cannot be guaranteed. In order to improve the drawback, PhaseCut [27,28] was presented by accurately separating the amplitude and phase variables on the basis of PhaseLift. Due to the lifting, the two methods result in the prohibitive computational costs of two-dimensional signal reconstruction. Recently, a convex relaxation PhaseMax was proposed [29][30][31][32][33][34] to deal with a PR problem with non-lifting and operating in the original signal dimension. It is worth noting that convex relaxation methods require multiple sets of amplitude measurements in the observation domain to ensure convergence.
The other category of methods is usually based on minimization or maximization of some cost function or fitness function by global optimization procedures [35][36][37][38]. In [35], the weighted L2-norm regularize as a multiplicative constraint is introduced in the cost function, and the conjugate gradient method is employed to minimize the cost function. In [38], the gradient descent algorithm is utilized to minimize the cost function. Due to a lack of measured phase information, the overall cost function in these methods is non-linear, which usually fail to obtain the global optimal solution. Furthermore, these methods highly rely on the accuracy of the initial phase guess. If initial estimation is not reasonable, it is very difficult to obtain relatively accurate phase information. In [36], the method introduces a suitable relaxation in the cost function and decomposes the cost function maximization problem into a limited number of different convex programming problems. Nevertheless, the accuracy of this method is closely related to the number of unknowns. In [37], the cost function is minimized by the genetic algorithm. While the cost function contains two functions to accelerate the convergence of the algorithm, the algorithm still needs nearly 100,000 times to converge. More than this, the acquisition of two sets of amplitude information on different closed cylindrical surfaces will result in a higher cost than on planes.
While a number of phaseless antenna measurement approaches have been developed, there are very few papers concerning the phaseless characterization of quiet zone in CATR. In [38], the technique used to retrieve the phase of the quiet zone in CATR is based on a given initial phase range and two sets of amplitudes on the planes. The area of the plane is much larger than the aperture area of the main reflector; meanwhile, the least sampling interval needs to be λ/4. Even so, the measured surfaces still result in a truncation error inevitably. It should be noted that when there are many unknown parameters to be solved, it will converge to a false solution. In [39], we use the alternating projection method to recover the phase. However, four sets of amplitude data need to be sampled. Therefore, a large number of unknowns involved and the strict requirements for achievable measurement accuracy make the related PR problem of large electrically sized CATR still a very difficult task. To this end, we introduce a new alternating projection method to retrieve phase information of the quiet in CATR without relying on the accuracy of the initial phase guess.
In this paper, the proposed alternating projection method is based on the combination of GS and HIO algorithms with additional spatial phase perturbation (SPP) which allows to escape from the attraction region of one local optimum to that of another candidate solution. The novel algorithm is implemented in three parts. The first and second parts are a combination of the improved GS and HIO algorithms with the SPP and the third part is a pure improved HIO algorithm. This method takes advantage of the fact that the introduced SPP can allow simply skipping from one attraction area to another attraction area by simply changing the added interference phase. Meanwhile, the introduction of HIO in the third part of the algorithm is an attempt to prevent excessive interference when the solution is near the optimal solution because excessive interference may result in falling into the current local minimum solution of the problem.
The performance of this method has been pointed out by extensive numerical analysis which uses synthetic data to retain all the features of a CATR. The data is generated using MATLAB and commercial software GRASP-10.0 developed by TICRA. Referring to a Cassegrain Gregorian tri-reflection CATR with a spherical main reflector and two shaped sub-reflectors introduced in [5], problems related to the practical application and performance of this system have been pointed out.

Phase Retrieval as a Feasibility Problem
The image recovery problem in a Hilbert image space L is to estimate the original form x based on the measurements of the related physical image and some a priori information [40,41]. In the phase recovery problem of this article, the measurements are composed of two sets of amplitude data, which are the modulus s of the original signal x from the object domain and the modulus m of the Fourier transformx of x, namely |x| = s and |x| = m. The original signal x ∈ C is complex value. The measurement data s and m are reals and nonnegative. Generally, we can properly model discrete physical signals/images in Hilbert signals/images space L. Thus, a signal x in L is a square-integrable function, mapping a discrete variable to a complex number.
It is necessary to consider an important piece of information that the support of x is included in some set D ⊂ R N . Moreover, the complement of D is specified as D. The characteristic functions of D and D are 1 D and 1 D , respectively, and the values of them are defined as S and M are defined as the signal sets in L satisfying the constraints of the object domain constraint and the Fourier domain constraint, respectively. The object domain constraint restricts x to the set S = {y ∈ L | y · 1 D = 0 and y 0}. ( Meanwhile, the projection of a signal x ∈ L onto the object domain constraint S is given by Since the measured amplitude s of the object domain is known, x can be defined by M is a signal set that satisfies the Fourier domain constraint, whose Fourier modulus is consistent with the modulus m, namely whereŷ denotes the Fourier transform of y, and is defined bŷ It should be emphasized that M is closed but non-convex, whereas set S defined by Equation (2) is a closed convex set (as seen in [42]). The projection of a signal x ∈ L onto the Fourier magnitude constraint set M is to replace the amplitude of the Fourier transform of x with the known measured data m and inversely transform the result, and is expressed as Therefore, the PR problem in Hilbert space can be expressed in mathematical form. This is to find a signal x ∈ L that satisfies the above constraint, namely Through the mathematical formula, the problem of solving PR can be transformed into finding a suitable point in the intersection of two constraint sets. This kind of problem is called a feasibility problem in mathematical optimization. In addition, in this article we will only consider the phase recovery problem in the case of S ∩ M = ∅. Because the sampled signals are discrete signals with finite dimensions, the Hilbert space in this article can be equivalent to a Euclidean space whose dimension depends on the amount of data required to deal with the problem.

The Proposed SPP GS/HIO-HIO Algorithm
Let us consider the PR of the quiet zone in a Cartesian coordinate system Oxyz of a triple offset reflector CATR. The phase retrieval problem of the plane S 1 at z = z 0 with the transversal dimension a * a in quiet zone can be solved when the amplitudes of the two planes S 1 and S 2 are measured as shown in Figure 1. The field of plane S 2 is the Fourier transform of the field of plane S 1 and this process can be achieved by a metasurface lens as explained in [43]. When alternating projection methods solve the feasibility problem of phase recovery, the signal is alternately projected onto the Fourier constraint set M and the object constraint set S to update until the algorithm converges and outputs the retrieved phase. If the update of the signal always satisfies the invariable constraint, the algorithm will not converge or fall into a local minimum corresponding to the "false solutions" of the problem. In order to avoid the shortages, the SPP GS/HIO-HIO algorithm proposed in this paper adaptively introduces different constraints to update the signal according to the process of dealing with the feasibility problem.
The process of the proposed SPP GS/HIO-HIO algorithm is shown in Figure 2. It can be mainly divided into the following three parts. The first part mainly uses the SPP GS algorithm to preprocess the signal. When the initial input signal consisting of a known amplitude and a randomly generated phase is projected onto the object domain, the constraint of SPP GS is performed to update the signal with the count number k increased by one. The second part is mainly to update the signal cyclically for k s − 1 times. Meanwhile, the constraint of SPP HIO is used to update the signal when it is projected onto the object domain under the condition of mod(k.k s ) > k 0 . The updated constraint of SPP HIO is replaced by that of SPP GS if mod(k.k s ) < k 0 , and then returned to the first part with the updated signal as the input value. The addition of SPP in the first and second parts is to prevent the algorithm from failing to converge. However, extra interference effects in the presence of SPP will occur when the updated signal is close to the original target signal, which will also lead the algorithm to fall into a local minimum corresponding to the "false solutions" of the problem. Therefore, SPP will be removed on the condition of k > K in the third part. Moreover, the constraints for the signal updating in the object domain are those of the HIO algorithm. The retrieved phase of the updated signal will not successfully generate until the algorithm converges.
In order to prevent the occurrence of overlapping effects at the boundary, we additionally introduced the zero-padding method by adding equal zeros in both the x-direction and the y-direction, introduced in the object domain. As shown in Figure 3, the amplitudes of plane S 1 can be divided into two sets, set D and its complement D according to the zero region. The elements of set D are sampled amplitudes, while the elements of D are all padded with zeros. Furthermore, the sampled region in plane S 2 should cover the entire zero-padding part.
On the basis of the zero-padding method, we present in Figure 4 the flow chart of the signal updating with constraints of the SPP GS algorithm and the SPP HIO algorithm, respectively, when the signal is projected onto the object domain. The flow chart of the HIO algorithm is not given for the reason that it only lacks the addition of the SPP compared with that of the SPP HIO algorithm.
Equation (9) is where the signal f k (x, y, z 0 ) is projected onto the Fourier domain constraint set M to obtain f k (x, y, z 0 ). Then the signal f k (x, y, z 0 ) is projected onto the object domain constraint set S and updated to f k+1 (x, y, z 0 ). According to the program process, the object domain constraint for signal updating will be different. When k ≤ K , SPP (e −j∆θ k ) is added and then e −jθ k+1 = e −j(∆θ k +θ k ) . If mod(k, k s ) ≤ k 0 , f k+1 (x, y, z 0 ) is defined as If mod(k, k s ) > k 0 , f k+1 (x, y, z 0 ) is defined as When k > K , SPP is removed and then e −jθ k+1 = e −jθ k ; the signal updating is updated to As explained in [12], the phase ∆θ k (x, y) in SPP is expressed as: where α is a scaling factor that determines the initial size of the perturbation. In order to ensure that the introduced disturbance allows the algorithm to avoid converging to a local minimum, the initial value of perturbation cannot be selected as too small. The default value of α in this paper is set to π/2 unless otherwise specified. Moreover, rand (−1,1) (x, y, z 0 ) is a randomly generated matrix at the kth iteration with the same dimensions as the plane S 1 being processed. The elements in the matrix are uniformly distributed between −1 and 1. The factor β is a constant feedback parameter. It has been found that the value of β between 0.5 and 1 can produce good results [44] and the value of β is set to 1 in this paper.

An Outline of the Solution Steps
The problem of PR has been transformed into a feasibility problem in mathematical optimization in Section 2. The main purpose for the optimization is to find a suitable signal in the intersection of two constraint sets. The following will explain in detail how the alternating projection algorithms work and the type of role that they play in the problem. For each update, the signal is first projected onto the Fourier-domain constraint set M and then onto the object-domain constraint set S. Furthermore, set S is convex (and actually linear) and closed, M is closed but nonconvex. The Fourier domain constraint remains the same throughout the solution. However, the support constraints contained in the object domain constraint set vary according to k.
When k ≤ K and mod(k, k s ) ≤ k 0 , the updated signal can be obtained through [45] f k+1 (x, y, Then, the signal can be further defined as For a better understanding, the projection process of the SPP GS method is shown in Figure 5. When k ≤ K and mod(k, k s ) > k 0 , the support constraints for the signal updating will be changed as (16), which can be rewritten as (17).
β is set to 1 in (17). Now let R S = 1 + e −j∆θ k P S − I and R M = 1 + e −j∆θ k P M − I. Thus, the updated signal can be rewritten as (18).
We can further simplify the form of (18) to obtain According to the form of (19), the corresponding signal projection and update can be divided into three steps through [46], shown in Figure 5. Firstly, calculate the reflection r k+1/2 of f k in relation to set M, and then calculate the reflection r k+1 of r k+1/2 in relation to the set S. f k+1 is the midpoint of the segment between f k and r k+1 . Finally, the updated f k+1 can be obtained from (19). If k > K , the perturbation phase e −j∆θ k will be removed, and the support constraint of the object domain set S for the signal updating is modified as (20), which can be simplified as follows: Let R S = 2P S − I and R M = 2P M − I; f k+1 (x, y, z 0 ) can be expressed as The corresponding signal projection and update through (22) has been shown in Figure 6. Similarly, r k+1/2 (r k+1 ) is the reflection of f k (r k+1/2 ) in relation to the set S (set M). f k+1 is the midpoint of the segment between f k and r k+1 . In this paper, the value of k 0 , k s , and K are set to 1, 40, and 600, respectively. The number of updates k starts from 1, and the total number of updates K is set to 1200. The whole projection procedure of the proposed algorithm under these parameters is shown in Figure 6.

Numerical Verification of the Theoretical Results
Several numerical experiments have been performed to evaluate the effectiveness and stability of the algorithm in the proposed Cassegrain Gregorian tri-reflection CATR [5]. The CATR consists of a spherical main reflector with a circular aperture, two shaped subreflectors and a Gaussian feed horn as shown in Figure 7. The aperture diameter of the main reflector is a (a = 1 m), and the diameter of the cross section of the effective quiet zone is b (b = 0.7 m). For the sake of simplicity, set the z-axis coordinate of the plane of the quiet zone to be retrieved in this paper to z 0 (see in Figure 3). The measurement planes S 1 and S 2 are sampled at λ/2 by commercial software GRASP-10 and MATLAB simulation software, respectively. If we want to obtain the missing phase of the plane with a diameter of a, the sampling areas of S 1 and S 2 are a * a and 2a * 2a, respectively. In order to assess the stability of the proposed approach, the missing phases of the S 1 with the CATR working at 75 GHz, 138 GHz, and 220 GHz, respectively, were recovered. Additionally, under the same initial input and total number of iterations (K = 1200), the recovery performance of the proposed algorithm is compared with that of the other four algorithms: GS, HIO, GS/HIO [12], and hybrid GS/HIO [47] algorithms.
Before describing the details of PR, some explanations on performance indicators are needed. The discrete sum squared error (SSE) between the amplitude of the Fourier transform of the current iteration F k (u, v, z 0 ) and the Fourier amplitude data m(u, v, z 0 ), is employed to judge the convergence of the algorithms and presented to be Moreover, the comparisons of the cost function values with different algorithms are given with the tri-reflector CATR working at different frequency bands. Considering the performance of the quantitative evaluation algorithm, we also define the phase normalized root mean square error (NRMSE) as: The phase of the plane S 1 located at z 0 was retrieved by the proposed algorithm and the other four algorithms for comparison. It can be found from Figure 8 that the retrieved size is the aperture size with a circle domain of 1 m. In order to simplify the comparison with different algorithms, only the normalized phase retrieved values on the horizontal (y = 0) and vertical (x = 0) axes are drawn. It can be seen that the performance of the proposed SPP GS/HIO-HIO algorithm is the best. Whatever frequency the CATR works at, the phase retrieved by the proposed algorithm can always be highly consistent with the target phase obtained by commercial software GRASP-10. The phase jitter trend recovered by the GS/HIO algorithm is roughly the same as that of the target phase. The performances of the other three algorithms are not good, especially at f = 75 GHz.
It can be seen from Figure 9 that the SSE decreases as the iterative number increases and the SSE performance converges when k = 1200 at different frequencies. Note that the lowest SSE performance is achieved by the proposed algorithm and that that of other algorithms no longer decreases or decreases very slowly when k > 600. The SSE performance of the proposed algorithm can further decrease sharply in the range of 600 to 800 because the constraint of SPP GS/HIO is changed to that of HIO when k > 600 to prevent the proposed algorithm from converging to the current local optimum. It is worth noting that the convergence performance of the proposed algorithm at f = 75 GHz is better than that at f = 138/220 GHz because the wavelength is larger at f = 75 GHz, and then phase jitter in QZ is larger due to the edge diffraction. In addition, fewer unknowns need to be solved in lower frequency. These factors will make the algorithm easier concerning obtaining large error values in the solution process, which is conducive to the algorithm escaping from the current local optimal solution to find a more appropriate solution in the derivative operations. In order to further compare the difference between the phase recovered by the algorithms and the target phase, the 2-D target phase distribution and the recovered phase distribution at f = 220 GHz are shown in Figure 10. Additionally, the NRMSE performances of phases recovered by each algorithm at different frequencies are listed in Table 1. In Figure 10a, the distribution of the target phase in the quiet zone on the y-axis is approximately symmetric, while the phase distribution on the positive half axis on the xaxis has more severe fluctuations than that on the negative half axis. Compared with target phase distribution of the plane with a diameter of 1 m, the phase distributions retrieved by the GS, HIO and hybrid GS/HIO algorithms, respectively, lack the main phase jitter details (see Figure 10d-f). The retrieved phase distributions by the three algorithms do not have a significant jitter in the quiet zone with a diameter of 0.7 m, but the jitter sharply increases outside the quiet zone. The phase distribution obtained by the GS/HIO algorithm contains too many phase jitters, especially in the quiet zone, while these jitters do not exist in the target phase distribution (see Figure 10a,c). As expected, the detail of the phase distribution retrieved by the proposed SPP GS/HIO-HIO algorithm is almost identical to that of the target phase distribution (see Figure 10a,b). A comparison of the computational time of the five algorithms under the same conditions has been shown in Figure 11. It can be seen that the GS algorithm has the fastest convergence speed, followed by the proposed algorithm. However, the GS algorithm easily converges to the wrong solution. Then, the performance of the proposed algorithm is the best from a comprehensive perspective.
Based on the analysis results of the above sections, the performance of the proposed algorithm is the best. In addition, the results of the phase presented in Figures 8 and 10 are retrieved based on the amplitude obtained under ideal conditions. In order to evaluate the robustness of the proposed algorithm, the phase will be retrieved under the interference of different degrees of additive white Gaussian noise at f = 220 GHz. We also use the GS/HIO algorithm to retrieve the phase under the same conditions for comparison. In order to facilitate comparison, the normalized retrieved phase on the horizontal axis and vertical axis with a range of the quiet zone are plotted in Figure 12.
Additionally, the NRMSE values of the retrieved phase with a range of 1 m as a function of the signal-to-noise ratio (SNR) are depicted in Figure 13.
As can be seen in Figures 12 and 13, the performance of the proposed algorithm is significantly better than that of the GS/HIO algorithm. Whatever the value of SNR, the NRMSE performance of the phase recovered by the proposed algorithm is greater than that recovered by the proposed algorithm. Besides, the phase recovered by the GS/HIO algorithm is in good agreement with the target phase in the condition of SNR > 30 dB. In Figure 13, when the SNR is greater than 25 dB, the NRMSE value of the phase recovered by the proposed SPP GS/HIO-HIO algorithm decreases sharply.
Notably, the missing phase of the plane with diameter a retrieved by the algorithms above is based on the amplitudes of planes S 1 and S 2 . The required amplitude sampling range of plane S 1 in the spatial domain is a * a with the remaining part padded with zero, and that of plane S 2 in the Fourier domain needs to be 2a * 2a. The sampling area of S 2 is too large and brings a huge burden of sampling and computing. Therefore, it is necessary to consider the sampling area of the S 2 plane. The following shows the comparison between the retrieved phase and the target phase of the plane with diameter a (a = 1 m), when the sampling area of the plane S 2 is reduced in different proportions with the sampling size still equal to λ/2. The results retrieved by the GS/HIO algorithm and the proposed algorithm are shown in Figures 14 and 15, respectively.   Comparing Figure 14 with Figure 15, it can be further verified that the proposed algorithm performed better than the GS/HIO algorithm. Furthermore, when the sampling area of S 2 is reduced to the same size as the plane to be retrieved, namely a * a, the phase retrieved by the proposed algorithm can be completely consistent with the target phase. If the sampling area of S 2 is further reduced, there is a difference between the recovered phase value and the target phase value (see Figure 15). It can be seen from the figures that the relative difference on the entire coordinate axis increases as the sample region decreases. Nonetheless, the relative differences between the recovered phase jitter and the target phase jitter are small in the region of a < 0.7 m. Generally, we mainly focus on the relative difference of phase jitter in the core region of the quiet zone and the difference between the maximum and minimum of phase jitter is the main evaluation criteria to be considered. Therefore, the proposed algorithm can still basically ensure high PR accuracy, when the sampling area of S2 is reduced to 1/16 of the original size.

Conclusions
This paper proposes an improved alternating projection method to solve the characterization problem of CATR with phaseless measurements. This problem has been formulated as a feasibility problem, that is, to find a suitable signal in the intersection of two constraint sets. The proposed algorithm contains three different update constraints, which perform a series of hierarchical optimizations of the problem, and the corresponding solution optimization process has been given. Meanwhile, the phase retrieval performances of other published alternating projection algorithms are compared with the proposed algorithm under the same conditions. The numerical examples are given in the Cassegrain Gre-gorian tri-reflection CATR condition and the performances prove that the method can retrieve the missing phase in a reliable manner without frequency limitation. Moreover, the algorithm can still guarantee higher PR accuracy when the acquired amplitude data is greatly reduced.