Phaseless computational imaging with a radiating metasurface

Computational imaging modalities support a simplification of the active architectures required in an imaging system and these approaches have been validated across the electromagnetic spectrum. Recent implementations have utilized pseudo-orthogonal radiation patterns to illuminate an object of interest---notably, frequency-diverse metasurfaces have been exploited as fast and low-cost alternative to conventional coherent imaging systems. However, accurately measuring the complex-valued signals in the frequency domain can be burdensome, particularly for sub-centimeter wavelengths. Here, computational imaging is studied under the relaxed constraint of intensity-only measurements. A novel 3D imaging system is conceived based on 'phaseless' and compressed measurements, with benefits from recent advances in the field of phase retrieval. In this paper, the methodology associated with this novel principle is described, studied, and experimentally demonstrated in the microwave range. A comparison of the estimated images from both complex valued and phaseless measurements are presented, verifying the fidelity of phaseless computational imaging.


Introduction
Remote sensing exploits the reflection of radiated waves to localize, image, and nondestructively detect objects under study. In imaging applications, measured waveforms are usually back-propagated to the object space using techniques such as Kirchoff migration [1], and Stolt's F-K migration [2], giving the location and reflectivity of the target of interest. This process requires the measurement of fast-time varying signals, represented by complex-valued harmonics in the Fourier domain. These phase measurements become challenging in the terahertz, visible, and X-ray ranges, and require the implementation of complex setups, often based on interferometric techniques. Thus, to overcome these hardware limitations, many phase retrieval techniques have been proposed in the last decade to allow for the reconstruction of complex field distributions solely from intensity measurements. They take inspiration from the concept of holograms defined by Gabor [3], applying alternating projection algorithms introduced by Gerchberg and Saxton [4] and Fienup [5,6], where at each iteration a complex field source is tailored such that the absolute value of its projection in the target space matches the measured intensity. In the last decade, a resurgence of attention in the scientific community has focused on the development of efficient phase retrieval algorithms. Independent and almost simultaneous contributions by Papanicolaou et al. [7] and Candès et al. [8] proposed comparable optimization techniques, both based on semi-definite programming of a relaxed problem. The practical implementation proposed by Candès et al. is presented here as the foundation of the approach proposed in this article. This technique focuses on the reconstruction of three-dimensional objects from the measured intensity of coded diffraction patterns (Fig. 1). A reconfigurable phase plate encodes the fields diffracted from a molecule illuminated by a coherent beam, allowing the reconstruction of a 3D electron density map from multiple intensity measurements of modulated projections. The depicted X-ray imaging setup can be considered under the emerging framework of coherent computational imaging achieved on the physical layer [9,10,11]. The aim of this paper is thus to highlight these similarities and demonstrate the reconstruction of near-field 3D images from phaseless measurements taken with a computational system. To this end, the mathematical derivations associated with the phaseless measurement of coded diffraction patterns are presented, followed by a concise review of recent phase retrieval techniques. The diffracted patterns measured on a plane made of m pixels y y y ∈ R m are expressed as follows: y y y = |A A Ax x x| 2 (1) where x x x ∈ C n is an unknown object and A A A ∈ C m×n is a known transfer matrix in which each line a a a i , i = 1, ..., m stands for a coded diffraction. The wave diffracted from the object is thus coded by a random mask, giving an illumination pattern y y y (l) of the form [12]: y y y (l) = |F F FD D D (l) x x x| 2 (2) where F F F ∈ C n,n is a discrete Fourier transform matrix and D D D (l) ∈ C n,n is a diagonal matrix whose entries are the known random complex transmittance of the mask modulating the diffracted pattern. L random masks are used for the estimation of x x x, leading to m = nL measured points to match the initial formulation of Eq. 1. An optimization problem is formulated to find the rank-1 matrix X X X = x x xx x x † . Indeed, for each line y k : where a a a i a a a † i is a rank-1 matrix. X being the outer product of two vectors, this matrix must satisfy: X X X 0, rank(X X X) = 1, y i = Tr(a a a i a a a † i X X X) for i = 1, ..., m dropping the rank constraint and replacing it by a trace minimization, accounting for the sum of the singular values of X X X. The semidefinite program PhaseLift is hence defined as: subject to X X X 0, Tr(a a a k a a a † k X) = y k , k = 1, ..., m The convexity of this formulation makes it solvable with the help of optimization software such as CVX [13]. Among the numerous applications of this approach, the retrieval from phaseless far-field data of a microwave array has recently been demonstrated in numerical studies [14,15]. This represents a promising approach for the simplification of radiation pattern characterization and far field imaging systems. However, the dimensionality "lifting" imposed by this algorithm represents the main drawback, squaring the number of unknowns to create the rank-1 matrix X X X.
Alternatively, novel alternating projection algorithms represent an interesting approach for solving large quadratic systems. They demonstrate the efficiency of iterative phase retrieval techniques, such as the block-Kaczmarz method [16], derived from the algebraic reconstruction technique [17], and the Wirtinger flow [18]. Recently, the truncated Wirtinger flow was proposed by Chen and Candès [12]. It has been demonstrated that these methods always converge to a solution when satisfying support constraints in the spatial domain and appropriate frequency oversampling, in contrast to the most popular methods introduced by Gerchberg, Saxton, and Fienup that can stall in local minima [8]. The truncated Wirtinger flow is the solution adopted in this article for its demonstrated efficiency, although the other mentioned methods remain compatible with quadratic formulations equivalent to Eq. 1. This algorithm computes the following equations for each iteration: For each iteration t, the value of x is thus updated by this descent gradient-like algorithm where µ t is a step size that can be determined for example by a backtraining line search and ∇ l i is a descent direction. The algorithm is computed on the adaptive index set S t+1 as determined by Chen and Candès [19], satisfying for any i ∈ S t : These constraints improved the efficiency of the initial truncated Wirtinger flow by dropping some gradient components ∇ l i (x) bearing too much influence on the search direction. The efficiency of this algorithm is demonstrated in [19] considering both ideal and noisy measurements-showing that the computational effort required for solving a random quadratic system is on the same order of magnitude of that of a linear system of the same size.

Phase retrieval adapted to near field imaging using a radiating metasurface
We propose the adaptation of the phase retrieval framework to a microwave computational imaging system by replacing the coding phase plate presented in the depicted X-ray system (Fig. 1) by a radiating metasurface and implementing the truncated Wirtinger flow instead of the PhaseLift presented earlier. The codes made openly available by their authors, Chen and Candès, are adapted to this study. Several studies have demonstrated the possible simplification of imaging hardware exploiting the recent development of spatially resolving antennas [20]. In this context, systems radiating structured illumination patterns have been proposed, exploiting the inherent structure in the imaged targets to limit the hardware complexity [9,10,21]. To this end, frequency diverse [11,22] and dynamically reconfigurable radiating apertures [23] were studied, demonstrating the linear dependency between the target space and the measured signals through tailored sensing matrices. In this paper, the principle of computational imaging is thus adapted to the phase retrieval problem using a frequency diverse radiating structure, taking benefit of the hardware simplification allowed by both approaches to propose a new paradigm for imaging systems. This demonstration is proposed in the microwave range as a proof of concept to easily compare the reconstructions from complex valued and intensity measurements, paving the way for millimeter wave, terahertz, and photonic applications.
Here, we study the mathematical formalism and develop the conditions in which intensity measurement of compressed waveforms can be applied to retrieve the positions of targets in the near field. The experimental setup and the associated parameters are defined in Fig. ( In this setup, a frequency diverse structure similar to those introduced in [9, 22, 24, 25] is considered. The large modal diversity excited by these metasurfaces allows for the radiation of structured field patterns which vary with the driving frequency. The quality factor of the metasurface is optimized to avoid modal degeneration, allowing for the radiation of a large number of pseudo-orthogonal spatial modes sensing the target space. The expression of the measured frequency signal ρ(ν) on the radiating device's output port can be expressed as: where φ (r r , ν) stands for the near field distribution of the metasurface measured at the aperture locations r r for each frequency ν, g(r r , r, ν) represents the Green function modeling the propagation of field from the object space r and the metasurface's aperture, and f (r) corresponds to a field source that is localized with this computational system. This problem can be expressed using a matrix formalism by discretizing equation Eq. 10; we represented the resulting vectors and matrices in bold notation as r r r = [r i ] 1≤i≤n , r r r r = [r r j ] 1≤ j≤n r , and ν ν ν = [ν k ] 1≤k≤m . A sensing matrix H H H ∈ C m×n is defined by the product of the Green function matrix G G G n,n r (ν k ) and the cavity near-field response written in the vector form φ φ φ n r (ν k ) for each frequency ν k : The sensing matrix allows for a representation of the linear dependency between the measured frequency signal ρ ρ ρ ∈ C m and the object f f f ∈ C n , leading to the following formulation: Previous works demonstrated methods of reconstructing the image vector f f f under certain invertibility conditions of sensing matrix H H H, corresponding to a sufficient number of radiated orthogonal patterns interrogating the target space [9]. This work extends the frame of computational imaging by considering the measurement of the intensity of the compressed waveform ρ ρ ρ , described by the following quadratic equation: Similarly to the coded diffraction problem, an object is reconstructed from phaseless data knowing the complex transfer function of the coding system. However, the frequency diversity exploited by this approach coupled to the associated derivation makes it distinguishable for its compatibility to near-field imaging, without using any actively reconfigurable radiating components. A simulation of the system depicted in Fig. 2 is studied to identify the number of independent modes required to ensure the reconstruction of the object under study.
In the numerical model of the frequency dependent and randomized radiation pattern of the metasurface, the radiating plane is set at y = 0, where y is the propagation axis. In accordance with (10), the radiation f (r r r) of a field source (set in the Fresnel zone of the radiating metasurface), is propagated to the aperture plane r r r r , compressed by its near-field responses φ (r r r r , ν ν ν) into a unique measurement at the port where the phaseless measurement is performed. Careful consideration of the dispersive nature of the metasurface φ (r r r r , ν ν ν) is required to study the convergence of the computational phase retrieval problem. In the numerical simulations, the impulse response is defined in the radiating plane by a Gaussian random distribution R d (r r r r ,t t t) with mean value zero, variance one, and discretized over steps of c/2ν max in space and 1/B in time. The quality factor Q of the metasurface is included in this model to consider the cavity's intrinsic losses and the coupling with all of the irises, leading to a modal degeneration [20]. The random field distribution is thus multiplied by a decaying envelope or damping time τ = Q/πν 0 , ν 0 being the central frequency of the studied bandwidth ( Fig. 3) [20, 21].
Then, computing the Fourier transform from the time domain to the discrete set of frequency components ν, the full model of the metasurface spatial and frequency response is: In coherent computational imaging system, where complex-valued signals are measured, dispersive antennas presenting long lasting responses are designed to limit the correlation between each radiation pattern in the frequency domain, improving the conditioning of the sensing matrix H H H. In this way, an estimation of the target signaturef f f can be computed followingf f f = H H H + ρ ρ ρ, where the superscript + stands for the Moore-Penrose pseudo-inverse operator [26]. An ideal metasurface would presented perfectly orthogonal radiation patterns ensuring that H H H + H H H = I I I.
In practice, metasurfaces are designed to have low correlation among patterns exploiting the frequency diversity, leading to a non-ideal inversion of the sensing matrix and the apparition of parasitic sidelobes. If we are limited to the measurement of intensity |ρ(ν ν ν)| 2 , such a direct approach can not be applied. But, alternating projection algorithms may be adapted to this problem to determine of the impact of a radiating metasurface's characteristics.
The spatial domain sampling is determined by the dimensions of the metasurface, modeled by an array of frequency dispersive dipoles with responses φ (r r r r , ν ν ν). According to the Rayleigh resolution limit, the transverse spatial sampling for a radiating array of dimensions D x and D z at a distance R is d x = λ c R/D x and d z = λ c R/D z . The range sampling for a wideband system is given by the width of the radiated pulses as d y = c/(2B), where B = ν max − ν min is the operating bandwidth.
Assuming the measurement of the intensity of a frequency signal ρ ρ ρ described by (13), a numerical study is carried out to estimate the criteria required for an accurate reconstruction of a spatial field distributionf f f I , where the subscript I denotes a reconstruction from an intensity measurement. The performance of such a computational system can be validated against a reconstruction of complex-valued measurements,f f f , with a relative error computed for each simulation as: Because the reconstruction from phaseless techniques is unable to estimate an absolute phase, a rotation of the samples by θ = f f f I ,f f f is performed to align the estimations in the complex plane before the subtraction. According to [16,18], successful estimations are considered when ε 10 −5 . The field source domain is defined by a number of voxels of positions r r r. The ratio between the number of measured modes m and the number of reconstructed voxels n, determines the size of the sensing matrix H(r r r, ν ν ν). The ratio m/n is considered in this study for different values of Q to determine the relation between the number of intensity measurements and the size of the reconstructed domain. The simulation is performed on a frequency band defined between ν min = 17.5 GHz and ν max = 26.5 GHz (K-band), with a frequency sampling dν = (ν max − ν min )/m.
The metasurface delay spread τ and the equivalent quality factor Q are defined according to the frequency sampling as: where α t is a frequency sampling parameter set according to the metasurface's damping factor τ. According to [20], the optimal sampling is α t = 1/π ≈ 0.32 considering that at most πτB orthogonal channels can coexist due to modal degeneration. The study is presented for α t ∈ [0.1, 2] to demonstrate the impact of frequency averaging on the phase retrieval algorithm. A domain of n = 6 3 = 216 voxels is considered in which a complex random field f f f must be retrieved. An array of 20 × 20 radiating dipoles spatially spaced at λ min /2 in both transverse directions is considered to simulate a metasurface whose frequency response is defined by (14) (Fig. 4). 100 trials with random metasurface responses and field distributions are computed for each couple of parameters to estimate an empirical set ensuring the accurate estimation offf f I (Fig. 5). According to this numerical analysis, the probability of successful recovery tends to reach its maximum when m/n 5 and α t 1 π (Fig. 6).   . Empirical success rate according to the sampling m/n. α t must be larger than 1/π to ensure that there is at least as many measured modes m as orthogonal channels available in the sensing matrix H to reconstruct n voxels.
As predicted in [19], a sampling of a least m = 5n ensures a good agreement between the estimations of the field source f f f with and without the phase information. Furthermore, the relation between the quality factor of the designed metasurface and the frequency vector ν ν ν must satisfy: Considering the case of a minimal sampling m = 5n and the definition of the frequency sampling dν = B/m, the maximum number of voxels imaged with this technique takes the following form: Having identified the critical designing and sampling parameters allowing for accurate estimation of the compressed, phaseless measurement, the impact of additive noise is studied in the next section to evaluate the performance of the algorithm in the context of near field computational imaging. The model given by Eq. 12 is modified to account for a Gaussian additive noise η η η: With a sampling m = 6n and for different values of signal-to-noise ratio (SNR), 100 trials are computed with random field distributions f f f and metasurface responses characterized by α t = 2 to study the convergence of the truncated Wirtinger flow adapted to this near field computational imaging problem (Fig. 7).  The value of α t is deliberately chosen to be large in order to speed up the numerical convergence. According to the sampling m/n, the algorithm converges on average in less than 200 iterations. For each value of SNR, the average µ and standard deviation σ d is computed (Fig. 8). According to the SNR value, the phase retrieval algorithm demonstrates its efficiency by leading to an average relative error µ of 1/ √ SNR. Finally, a convergence study is presented considering the impact of α t acting on the frequency oversampling. A set of numerical simulations are computed with SNR = 10 6 and sampling m = 6n, considering 1000 trials for each value of α t (Fig. 9). This numerical study depicts the relation between the number of iterations required to reach convergence with respect to α t , demonstrating the positive impact of a high quality factor and a fine frequency sampling on the convergence speed.

Pratical implementation and experimental results
The theory of phaseless computational imaging is experimentally validated using a metasurface operating in the microwave regime. To this end, a metallic leaky cavity of 28.5 × 28.5 × 15.2 cm 3 was created (Fig. 10), inspired by a computational imaging prototype introduced in [22]. The front plate is perforated by a 15 × 15 cm 2 square array of circular holes randomly set on a 0.6 cm uniform grid. An open-ended waveguide source is set in front of the cavity and localized using the computational system. In contrast with the setup depicted in Fig. 2, two ports are connected to the back of the cavity (Fig. 11). In this way, different superpositions of modes can be measured by the two ports according to their locations in the cavity, increasing the amount of independent information in the frequency domain. A spherical reflector is also set in the cavity, providing additional mode mixing and increasing the number of uncorrelated states in the cavity [22,27].  To match the simulations, the operating frequency range is defined in the K-band between ν min = 17.5 GHz and ν max = 26.5 GHz, sampled by m ν = 3601 frequency points. In this way, m = 2m ν points are measured for the phaseless estimation of the n voxels constituting f f f . The patterns radiated at each frequency and for each excitation port are measured by a near-field, single-polarized probe moved on a planar synthetic aperture by a translation stage. This field is numerically back-propagated to the metasurface plane to estimate Φ 1 (r r r r , ν ν ν) and Φ 2 (r r r r , ν ν ν), the transfer functions of the cavity when exciting ports 1 and 2, respectively. Examples of near-field distributions are presented in Fig. (12) for two consecutive frequencies of the vector ν ν ν. Different pseudo-random spatial fields distributions are thus obtained as a function of the excited port and of the frequency, due to the low level of loss of this cavity. Fig. 12. Comparison of the near-field distributions Φ 1 (r r , ν) and Φ 2 (r r , ν) measured for the independent excitation of ports 1 and 2. The results are depicted for two consecutive frequency ν 1 = 23 GHz and ν 2 = 23.002 GHz of the frequency vector ν ν ν.
The measured quality factor of the cavity is about 12000 for both ports, determined by fitting exponential functions to the measured radiation patterns expressed in the time domain using a Fourier transform and taking the root mean square over the spatial dimension r r r r . A formulation matching Eq. 1 is proposed according to the measured signals ρ ρ ρ 1 ∈ C m ν and ρ ρ ρ 2 ∈ C m ν on ports 1 and 2 by concatenating the measured signals and the corresponding sensing matrices computed with Eq. 11: where ρ ρ ρ ∈ C m and H H H ∈ C m×n , dimensioned such that the target field distribution f f f ∈ C n is estimated, with m = 2m ν . The upper bound of the number of reconstructed voxels can be computed according to the quality factor and the setup presented in this experiment. Merging Eqs. 16 and 17 gives: From the expression of the time constant τ, the maximum number of orthogonal channels is bounded by [20]: In the considered experiment, there are at most m = 2m ν = 9800 independent frequency points, measured on the two ports of the cavity. Under the assumption of an ideal estimation of the sensing matrix H H H and a sampling n = m/5, a maximum of 1960 voxels can be reconstructed. For this validation, since m ν = 3601 frequency points are measured, a maximum of n = 1440 voxels can be estimated with a sampling of m = 5n. Considering this setup, the value of the parameter α t is given for one port by Eq. 17: Under these conditions and considering that two signals are measured for the estimation of f f f I , a fast convergence of the iterative phase retrieval algorithm is expected.
A first validation is proposed for a scene made of 10×10×10 = 1000 voxels centered around x = 0, y = 0.4 m, z = 0. In accordance with the numerical study, the spatial sampling is defined as d x = λ c R/D x ≈ 3.6 cm, d z = λ c R/D z ≈ 3.6 cm and d y = c/(2B) ≈ 1.7 cm, with R = 0.4 m, D x = D z = 0.15 cm, and B = 9 GHz. A probe is first set in front of the radiating metasurface in the middle of the pre-defined region of interest. A comparison between the retrieved field distributionsf f f andff f I is presented in Fig. 13.
The field distribution reconstructed from intensity-only measurements matches the estimation from the complex measurements, taken as a reference in this study. A higher level of noise is observed in the phaseless case, most likely due to a non-ideal estimation of the sensing matrix H H H considering that a supplementary mounting structure (source of diffraction) was added after the near-field characterization of the leaky cavity leading to a relative error ε = 0.57. A more precise comparison of the two estimated fieldsf f f andff f I is proposed, extracting the x, y, and z-cuts from the maximum values (Fig. 14).  Comparable locations and resolutions are observed in both cases, validating the fidelity of the truncated Wirtinger flow applied to this phaseless computational system. A larger domain is considered for the last part of this experimental validation, studying the impact of the sampling m/n presented earlier in a practical situation. A domain of 20 × 20 × 10 = 4000 voxels is thus considered this time, conserving the same spatial sampling and centered at the same location. We now consider a sampling of m/n = 9800/4000 = 2.45 (with the approximation of two independent measurements ρ ρ ρ 1 and ρ ρ ρ 2 ). In contrast with the numerical simulations where spatial random field distributions were considered, the experimental cases are focused on the reconstruction of a punctual point like object. Even if a relative error of ε < 10 −5 may not be reachable, we are interested in determining whether a localization of the field source is possible in the given conditions. A comparison of the reconstructions achieved with and without the phase information is once again presented from the same measurements as in the previous case (Fig. 15). A relative error of ε = 0.82 is obtained for this phaseless reconstruction, which is consistent with expectations of a larger ε than in the previous case due to the lower sampling m/n in this case. Despite this larger error, the localization of the field source remains possible, as depicted in the comparison between the field cuts presented in Fig. 16. Taking benefit of this larger domain, the source is translated at an off-center location to be imaged (Fig. 17). With a relative error of ε = 1.03, the three-dimensional reconstructions with and without the phase information remain comparable. The x, y, and z-cuts are extracted from the maximum value of both reconstructions for a finer analysis (Fig. 18).  While the fields extracted from the transverse axis are almost identical, a discrepancy is noted in the range axis. Indeed, the range information is coded by the time of arrival of propagating waves, mainly represented by phase ramps in the frequency domain that cannot be exploited in intensity measurements. Despite the lack of phase data and the under-sampling of the considered case, it remains possible to estimate the location of the field source with a degraded resolution compared to the computational imaging case based on complex valued measurements.

Conclusion
An application of a phase retrieval algorithm to a computational imaging system has been presented, allowing for the spatial reconstruction of field distributions from phaseless measurements. The truncated Wirtinger flow has been adapted in this study to determine the position of field sources from the measurements of a metasurface radiating pseudo-orthogonal patterns in the frequency domain. In contrast with the coded X-ray diffraction experiments simulated by Candès et al., there is no need of a reconfigurable random lens since the information is coded in the frequency domain by a static and passive device. While this application has been presented in the microwave range to where it is possible to compare the results with and without the phase information, the most useful applications stand at higher frequencies where the burdensome measurement of phase is problematic and limits the realization of 3D imaging systems. Future studies will thus focus on extending the implementation of such a technique to the terahertz, visible and infrared domains.

Aknowledgements
This work was supported by the Air Force Office of Scientific Research (AFOSR, Grant No. FA9550-12-1-0491).