## Abstract

Passive microwave imaging of incoherent sources is often approached in a lensless configuration through array-based interferometric processing. We present an alternative route in the form of a coded aperture realized using a dynamic metasurface. We demonstrate that this device can achieve an estimate of the spectral source distribution from a series of single-port spectral magnitude measurements and complex characterization of the modulation patterns. The image estimation problem is formulated in this case as compressive inversion of a set of standard linear matrix equations. In addition, we demonstrate that a dispersive metasurface design can achieve spectral encoding directly, offering the potential for spectral imaging from frequency-integrated, multiplexed measurements. The microwave dynamic metasurface aperture as an encoding structure is shown to comprise a substantially simplified hardware architecture than that employed in common passive microwave imaging systems. Our proposed technique can facilitate large scale microwave imaging applications that exploit pervasive ambient sources, while similar principles can readily be applied at terahertz, infrared, and optical frequencies.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Traditional imaging systems use lenses to achieve a direct mapping between an object and its image by physically manipulating fields of arbitrary coherence. In such systems, the effects of optical or electromagnetic field coherence can significantly impact the resulting image [1], manifesting as speckle or restricted depth of focus, for example. Nonetheless, the physical image-forming operation remains largely unchanged for different degrees of coherence. In contrast, lensless imaging and the associated computational approaches can be extremely sensitive to field coherence properties. Holographic image formation typically demands active illumination with fields of sufficient coherence to realize measureable interference effects. Similar requirements hold for methods such as computational ghost imaging [2–5] and array-based methods [6] that leverage complex electric field relationships for image formation through backpropagation and related techniques.

In the cases of partially coherent or incoherent fields, the relevant quantity is the coherence function, which characterizes the field’s spatiotemporal second-order correlations [7,8]. The coherence function is known to obey the wave equation for fields of arbitrary coherence [8], and so is amenable to computational manipulation through standard methods. In this way, passive characterization and imaging of primary or secondary sources can be achieved through remote measurement and free-space backpropagation of the coherence function [9–12].

In optical coherence imaging, one estimates an image of the source through samples of the coherence function measured using, for instance, an astigmatic coherence sensor [13]. Phase space tomography retrieves images from projections of a transformed representation of the coherence function [14,15]. Imaging correlography [16] and memory effect imaging [17,18] involve measuring far-field intensity correlations of coherent fields followed by a phase retrieval procedure and a Fourier transform to obtain the source coherence function. Above all, the most straightforward approach to passive characterization of a source distribution is through direct interferometric techniques, as epitomized by the Michelson stellar interferometer [7] and more broadly applied using microwave and RF technologies in radio astronomy applications [19].

Long-wavelength microwave and RF radiation is particularly robust to scattering by atmosphere and nonmetallic materials, making it well suited for applications such as nondestructive testing [20], biomedical monitoring [21], threat detection [22,23], satellite-based earth observation [24], and ground-based radio astronomy [19]. This property paired with mature measurement hardware enables accurate interferometric measurements by direct electronic signal mixing. The most successful application exemplifying this technique is array-based radio astronomy, depicted schematically in Fig. 1(a). In a typical radio astronomy application, one wishes to reconstruct the intensity distribution of an incoherent, quasimonochromatic source located in the array far-field. Under these assumptions, the source coherence function propagates to the antenna array according to the Van Cittert–Zernike theorem [7], which gives a Fourier transform relationship between the source intensity distribution and the coherence function measured at the array. A single element of the coherence function is measured by simultaneously collecting radiation at two distinct antennas and mixing the acquired signals physically, using a correlating circuit, or computationally through digital processing. Under the assumption of ergodicity, the correlation is evaluated with respect to a time average to achieve a sample of the coherence function. The diffraction-limited estimate of the coherence function is achieved by repeating this correlation procedure for each pair of array elements. The Van Cittert–Zernike theorem then gives a mapping between each cross correlation and an element of the source’s Fourier transform, allowing an image to be retrieved through an inverse Fourier transform of the measured coherence function.

In the approach just described, adequate interferometric characterization of the coherence function demands performing correlations between all pairs of array elements, which can result in unwieldy hardware and calibration, in addition to lengthy measurement times [23,25–27]. For some applications, these requirements may be alleviated through compressive imaging strategies [15]. Additionally, aperture synthesis techniques can substantially reduce the required number of array elements and associated hardware complexity [19]. In this work, we propose further reduction of the number of signal chains through a procedure fully analogous to the computational methods pioneered in coded aperture imaging schemes [28]. The coded aperture passive imaging concept is illustrated in Fig. 1(b). In its simplest form, coded aperture imaging is achieved by modulating an image with a series of spatially diverse patterns, and measuring the resulting intensities at as low as a single pixel [29]. In single-pixel imaging, each measurement corresponds to the projection of the image onto a given modulation pattern, so that an image can be formed as a linear superposition of the patterns weighted by the respective measurements [29]. More advanced reconstruction techniques involve inversion of the linear matrix equation describing the system forward model achieving, for example, compressive single-pixel [30] and spectral imaging [31]. In the following, we will show that one can apply similar principles in order to estimate a source distribution from samples of the coherence function multiplexed according to a linear matrix equation defined in terms of complex-valued modulation patterns and single-pixel measurements. In this case, the cross correlation operation reduces to the autocorrelation of a signal measured at a single position. These values can equivalently be obtained from power measurements in the frequency domain, a fact that we leverage to minimize phase sensitivity and hardware complexity in our system.

Spatial modulation necessary for aperture encoding can be achieved at microwave frequencies using electrically large chaotic cavities or metasurface apertures [32,33]. Metasurface apertures consist of a metallic waveguide structure with subwavelength, resonant metamaterial elements etched into the radiating surface. Operating in transmission or reception roles, these metasurface apertures achieve diverse gain patterns either by leveraging the dispersive response of the metamaterial elements to frequency variations in the guided wave [34], or by directly modifying the response of each element through electronic biasing [35]. The latter approach defines a dynamic metasurface aperture (DMA). In this way, metasurface apertures can multiplex spatial and spectral information onto a series of sequential single-port measurements.

Using frequency-diverse chaotic cavities at microwave frequencies, several works have demonstrated passive imaging of noise and thermal sources through computational principles under a radiometric formalism. In Ref. [36], the authors perform passive synthetic aperture interferometric radiometry (SAIR) through compressive encoding of a noise signal onto only two receiver ports using a mode mixing metallic cavity. The work in Ref. [37] builds on these efforts, outlining a matrix formulation for the decoding and near-field backpropagation steps. Computational passive imaging of true thermal sources is demonstrated in Ref. [38] using a similar chaotic cavity structure. While all of these works achieve passive imaging of radiating structures, the static nature of the encoding devices limits the accuracy of the reconstruction. To remedy this, the authors leverage frequency diversity of the cavities to improve the image intensity estimates, which achieves satisfactory imaging at the cost of spectral information. In addition, the radiometric formulation employed in these works emphasizes thermal sources over general noise sources, which can be treated under a coherence theoretical framework. The relationship between radiometric and coherence theoretical formalisms is well known [39,40], and provides a thorough framework for passive imaging of RF sources of incoherent radiation.

In this work, we demonstrate passive near-field spectral imaging by encoding the microwave coherence function using two-dimensional DMAs. We illustrate how dynamic modulation of the DMA encodes source spatial and spectral information onto single-port measurements. We further demonstrate the retrieval of spectral images using two acquisition schemes: spectral filtering by the acquisition hardware, and spectral multiplexing by the metasurface itself [41,42]. The DMA and its functionality are described in Section 2. In Section 3, we formulate the forward and inverse problems under a coherence theoretical framework, and experimentally demonstrate single-port, passive imaging of unknown sources in Section 4. We examine the effect of the encoding/decoding procedure as well as the near-field point spread function in Section 5 and conclude with discussion and suggestions for future studies.

## 2. DYNAMIC METASURFACE APERTURES

A waveguide-fed, DMA was first reported in Ref. [35] in the form of an effectively one-dimensional microstrip waveguide antenna. Recently, a two-dimensional printed cavity dynamic metasurface was introduced [43], which we employ for our studies in this work. Design and fabrication details are given in Ref. [44], but we repeat the information relevant to our purposes here.

The DMA used in this work is a printed cavity fabricated using standard
printed circuit board (PCB) technology (Figs. 2 and 3). A via
fence and two copper layers form the boundary of the cavity volume
composed of a 60-mil-thick Rogers 4003 dielectric substrate, and an RF
coaxial connector near the center of the cavity samples the diverse cavity
modes. The upper conductive layer is fashioned with 96 complementary
electric-LC (cELC) metamaterial elements [45] with Lorentzian resonant frequencies between 20–24 GHz. The
metamaterial radiators are divided into three distinct groups of 32
elements, each group defined by a different resonance frequency. The three
groups are distributed in a pseudorandom configuration over the whole $14 \times
14\;{\text{cm}^2}$ aperture area. The elements with the
highest resonance are approximately 3.5 mm in size. In addition, MACOM
MADP000907-14020W PIN diodes are incorporated into each element and
controlled through an external bias voltage, allowing independent on/off
control of each element. Individual elements are addressed through vias
connecting each element to a DC bias circuit constructed on the back of
the cavity layer. A DC bias signal is addressed to each element
independently using a series of daisy-chained shift registers controlled
by an Arduino microcontroller. A pair of radial stubs are used to decouple
the RF signal from the DC bias signal. We refer to the $m$th distinct bias voltage sequence as a
*tuning state*. These metamaterial elements
behave as magnetic dipoles described by a polarizability ${\alpha _m}({\vec
r_i},\omega)$ [46–48], which form discrete dipole moments ${\vec m_m}({\vec
r_i},\omega)$ by coupling to the $z$ component of the cavity magnetic mode ${H_z}({\vec
r_i},\omega)$ probed by the single RF port,

*modulation function*. In this work, we characterize the modulation function from scattering parameter measurements taken with a vector network analyzer (VNA), obtained by scanning in the radiative near-field with an open-ended waveguide (OEWG) probe polarized in the $\hat y$ direction [49,50]. The plane over which these patterns are characterized defines the aperture plane throughout this work.

In the following, we employ simple random tuning to retrieve an estimate of the source distribution. Optimal encoding of both spatial and spectral information by the resulting speckle modulation patterns would satisfy the condition ${\langle A_m^*(\vec r_a^\prime ,{\omega ^\prime}){A_m}(\vec r_a^{{\prime \prime}},{\omega ^{{\prime \prime}}})\rangle _m} = \delta (\vec r_a^\prime - \vec r_a^{{\prime \prime}})\delta ({\omega ^\prime} - {\omega ^{{\prime \prime}}})$, where $\langle\cdot\rangle_m$ denotes an average over $m$ tuning states. In practice, this correlation has a finite width among other artifacts resulting from aperture design, as described in Section 5. Using these modulation functions, the encoding of the $y$ component of the electric field incident on the aperture, ${E_a}(\vec r_a^\prime ,\omega)$, onto a single-port measurement ${g_m}(\omega)$ can be described mathematically as

Figures 2(b)–2(d) depict three experimentally measured magnitudes of the DMA modulation function, corresponding to the $y$ component of the radiated electric field, and reveal the receive pattern diversity achievable both through frequency variation and dynamic tuning.

The metamaterial elements are distributed sparsely over the DMA surface with an average spacing of about one free-space wavelength. The coupling of these elements to a single RF port thus allows the metasurface to perform effectively as a sparse 96-element dipole array. For comparison, at a frequency of 20 GHz, a dense array sampled at half a wavelength would consist of approximately 350 separate antennas. Neglecting aperture synthesis techniques, correlating the signals from every antenna pair in such an array would require 61,075 measurements. In contrast, we will demonstrate that our approach can achieve imaging through compressive reconstruction from only 1000 coded DMA measurements.

## 3. FORMULATION

In this section, we develop a model for the interaction of a target source distribution with our receiving DMA system. This model can describe microwave sources, such as antennas, that radiate unknown signals within our operating bandwidth. Given sufficient signal strength or integration time, sources of thermal radiation are also contained within our model, as are objects which scatter ambient signals from uncooperative sources and thus behave as secondary sources. Regardless of the exact physical nature of the source fields, we can characterize such targets as random, incoherent distributions from the perspective of our receiving DMA system, as long as these source fields are suitably uncorrelated both temporally and spatially. Our goal here is to formulate a method for passive image reconstruction of general sources obeying these assumptions. After describing our model, we detail the reconstruction algorithm, which includes a compressive decoding and backpropagation procedure.

We consider a random source distribution that is spatially and temporally incoherent. A single linear polarization component of the fields resulting from this source can be characterized by its scalar spectral density ${S_s}({\vec r^\prime},\omega)$. According to [1,8], this spectral density can be defined as an expected value over an ensemble of monochromatic field realizations $\{{E_s}(\vec r,\omega)\}$ for each frequency,

#### A. Independent Measurements

The aperture fields incident on the receiving DMA are modulated according to the complex-valued receiver modulation function ${A_m}(\vec r_a^\prime ,\omega)$. A single frequency component of the signal ${g_m}(\omega)$ measured on the DMA port for tuning state $m$ is then given by the inner product of Eq. (3). Using this relation, the average spectral magnitude measured at the single port of our receiving DMA, neglecting measurement noise, is given by

We find that the expected values of a set of spectral measurements are linearly related to the aperture CSD ${W_a}(\vec r_a^\prime ,\vec r_a^{{\prime \prime}},\omega)$ through Eq. (8). Namely, each measurement is the result of weighting the CSD by the respective modulation functions and summing over all pairs of positions.

Under our assumption of an incoherent source, we may express the full forward model as

In discrete form, the forward model for each frequency component can be expressed as

Here, ${{\bf g}_\omega}$ is an $M \times 1$ vector containing the spectral magnitude measurements, ${[{{\bf g}_\omega}]_i} = \langle |{g_i}(\omega {)|^2}\rangle$, and ${{\bf s}_\omega}$ is an $N \times 1$ vector that represents the spectral density distribution discretized into $N$ voxels, i.e., ${[{{\bf s}_\omega}]_j} = {S_s}(\vec r_j^\prime ,\omega)$. ${{\bf H}_\omega}$ is an $M \times N$ matrix with elements given by ${[{{\bf H}_\omega}]_\textit{ij}} = {H_i}(\vec r_j^\prime ,\omega)$.

In this form, estimating the source spectral density becomes a problem of inverting the ill-conditioned linear matrix Eq. (11), which can be achieved through a variety of well-known methods [37,51]. In this work, we assume a sparse source distribution and achieve compressive reconstruction as [52]

#### B. Multiplexed Measurements

The measurement model given in the previous section is appropriate when spectral measurements can be acquired using, for example, a spectrum analyzer with frequency-selective filters or time-domain waveform sampling. Meanwhile, the metasurface design employed here exhibits frequency diversity [43,54,55] in addition to the pattern diversity that results from dynamic modulation. If the different frequency components are encoded with sufficient independence, we can expect to be able to retrieve spectral information from multiplexed, frequency-integrated measurements [41,42]. This approach is analogous to retrieval of a temporal signal from time-integrated measurements in temporal ghost imaging [56,57].

The total intensity measured at the DMA port is given by the frequency-integrated expression

This can similarly be written in discrete form as

where ${[{\bf g}]_i} = \int_\omega \langle |{g_i}(\omega {)|^2}\rangle \text{d}\omega$ consists of the $M$ frequency-integrated intensity measurements, and the matrix ${\bf H}$ is formed by concatenating the matrices ${{\bf H}_\omega}$, corresponding to each of the ${N_f}$ frequencies, to form an $M \times (N*{N_f})$ matrix. ${\bf s}$ is constructed similarly to give an $(N*{N_f}) \times 1$ vector. Given a set of $M$ frequency-integrated measurements, we can obtain an estimate of the discretized source spectral density through compressive inversion of Eq. (15).## 4. EXPERIMENTAL RESULTS

As discussed in Section 3, source reconstruction can be achieved from a series of measurements multiplexed onto single-port magnitude measurements, a procedure which is largely insensitive to experimental factors such as cable length or measurement timing. We use this fact to build a simple setup for experimental verification of our proposed imaging system in the following. Our experimental setup is shown in Fig. 3, and consists of a single receiving DMA connected to an Agilent Technologies N9344C spectrum analyzer. The DMA is modulated using an Arduino microcontroller, and spectral magnitude measurements are taken over a frequency range of 18–20 GHz for each of 1000 tuning states by the spectrum analyzer operating in swept mode. To achieve an average spectral measurement as demanded by Eq. (8), the integration time was selected to sufficiently exceed the inverse of the filter bandwidth at each frequency point. Measurements are automated and coordinated using MATLAB. The total measurement time was approximately 45 min, and was determined largely by the signal integration time. Acquisition time can potentially be improved through the use of centralized, dedicated hardware. The sources to be imaged consist of various radiating structures emitting swept-frequency signals generated by network analyzers (Agilent Technologies N5245A PNA and E8364C PNA). The signals feeding the structures are in no way coordinated with the receiving hardware, and so behave essentially as noise sources radiating over the measured bandwidth. This is analogous to practical scenarios in which one wishes to passively image RF sources arbitrarily fed by uncooperative radio hardware.

We first demonstrate the ability to achieve passive imaging of a source intensity distribution by integrating the estimated spectral information. Imaging performance of extended incoherent objects is illustrated in Figs. 4 and 5. The sources were fed by approximately 5 dBm swept-frequency signals over the range 18–20 GHz. Measurements were obtained as described above, using a 3 s integration time to measure the entire bandwidth. Preprocessing consisted of binning the spectral measurements into 11 frequency intervals and summing the measurements within each interval in order to achieve the desired coarse frequency sampling. We then subtracted the mean over all DMA tuning states from each measurement, which served to center the data and associated measurement noise, and was empirically observed to improve contrast. An image was reconstructed for the mean frequency of each of the 11 intervals according to Eq. (12). The complex near-field scans comprising the modulation function ${A_m}(\vec r_a^\prime ,\omega)$ were additionally mean-subtracted before computing each matrix ${{\bf H}_\omega}$, which improved conditioning by minimizing correlations induced by coupling to the DMA coaxial connection. Finally, the total intensity of the source was computed by summing over all frequencies according to Eq. (13).

In Fig. 4, the source consists of the frequency-diverse metasurface aperture detailed in Ref. [58] and shown in Fig. 4(a), which is designed to radiate with pseudorandom radiation patterns over a frequency range of 17.5–26.5 GHz. The radiating elements on this metasurface are arranged predominantly along two stripe regions, as is evident in the image of Fig. 4(b) recovered using our compressive strategy to estimate the source distribution at a distance of 40 cm over the operating bandwidth of 18–20 GHz. For comparison, Fig. 4(c) illustrates a raster scan of the radiating metasurface target obtained from coherent measurements with a VNA. The scan image was convolved with the system point spread function evaluated at the mean frequency of 19 GHz using methods described in Section 5 to emulate the effect of diffraction limitations. In both cases, the radiation from individual slots is not distinguishable due to the finite resolving capabilities of the imaging system. As a result, the image consists of two strips, reflecting the arrangement of slots on the metasurface.

Figure 5 illustrates coded aperture imaging of a radiating cavity, in which the radiating surface has been perforated with slots in the pattern of a “smiley” face [Fig. 5(a)]. The radiating source is approximately 40 cm from the DMA receiver. Figure 5(b) depicts the reconstruction achieved by directly solving for the total source intensity using combined frequency information over the bandwidth 18–20 GHz.

Spectral imaging capabilities are demonstrated experimentally in Fig. 6. In this case, spectral filtering is performed by the spectrum analyzer, allowing an image to be formed from each measured frequency component. The sources consist of two OEWGs fed by separate VNAs that are not coordinated with each other. The OEWGs are located approximately 30 cm from the DMA aperture, and are offset longitudinally from each other by 5 cm. The first source is positioned at about $({y_1},{z_1}) = (- 0.05\;\text{m}, - 0.05\;\text{m})$, and radiates with a $-5 \; \text{dBm}$, swept-frequency waveform over the range 19–20 GHz. The second source is positioned at approximately $({y_2},{z_2}) = (0.05\;\text{m},0.05\;\text{m})$ and radiates with a 5 dBm, swept-frequency waveform over the bandwidth 18–19 GHz. Reconstruction consisted of preprocessing similar to that employed above, except that in this case the minimum over all masks was subtracted to remove DC offset while maintaining accurate spectral information. The representation in Fig. 6(a) was achieved by summing the spectral image over the first half of the operating bandwidth, resulting in an image of the blue point alone, while the orange point was similarly reconstructed by summing over the latter half of the operating bandwidth. Figure 6(b) gives a more quantitative representation of the estimated spectra. Here, the blue and orange curves correspond to the estimated spectra for the positions of maximum intensity in the blue- and orange-colored images of Fig. 6(a), respectively. Each spectrum was normalized to its maximum value prior to plotting. Deviation of the estimated spectra from the true radiated spectra is a result of the imperfect estimate achieved by our solution. The application of more sophisticated inversion methods can potentially improve this result, while the qualitative agreement obtained here indicates that accurate spectral information can be recovered over full, three-dimensional regions using the current method.

Next, we illustrate the ability to form spectral images from frequency-integrated intensity measurements multiplexed by our dispersive DMA. The multiplexed data was obtained in this case by summing the spectral magnitude measurements for each tuning state over the entire operating bandwidth. Spectral images were then obtained through compressive inversion of the matrix Eq. (15).

Figure 7(a) shows the result of spectral image reconstruction of an OEWG source from multiplexed measurements. We see that we are able to successfully localize this source at different transverse locations, while simultaneously retrieving an estimate of its radiation spectrum, as shown in Fig. 7(b). For this result, the measured spectra were truncated to different spectral profiles in postprocessing before performing the reconstruction. We compare the spectral estimates obtained from frequency-multiplexed measurements to those obtained by reconstruction with independent frequency measurements. The recovered spectral profiles largely agree, besides errors introduced by frequency correlations in the DMA response.

Finally, we demonstrate the ability to distinguish the spectra of simultaneously radiating sources using frequency-multiplexed data in Fig. 8. The sources in this case are the same as those reconstructed in Fig. 6. Figure 8(a) depicts the result of summing the recovered source distribution over the frequency range 18–19 GHz, which isolates the point located at approximately $(0.05\;\text{m},0.05\;\text{m})$. The other source at $(- 0.05\;\text{m}, - 0.05\;\text{m})$ is isolated and imaged by summing over the frequency range 19–20 GHz. The recovered spectra are shown in Fig. 8(c), and are seen to closely correspond to those estimated using independent frequency measurements as in Fig. 6.

## 5. ANALYSIS

In this section, we first detail some effects of the encoding procedure on the retrieved spectral density estimate. We then examine near-field imaging characteristics for our planar acquisition geometry. The system model described by Eqs. (15) and (10) indicates that the success of our imaging strategy and the conditioning of the measurement matrix ${\bf H}$ will depend largely on the degree of independence between our collection of modulation patterns. Figure 9(a) depicts the spatial autocorrelation of our mean-subtracted modulation function ${A_m}(\vec r_a^\prime ,\omega)$, computed over 1000 distinct tuning states at a frequency of 19 GHz for the center position of our aperture. An ideal encoding function would approximate a delta function, though the nonzero width of the autocorrelation in Fig. (9) reflects the finite resolving capabilities of our aperture. The sparse coverage of the DMA aperture with radiating metamaterial elements is seen to manifest as sidelobes in the autocorrelation, which leads to pseudorandom signal mixing in the spectral density estimate. Figure 9(b) illustrates the low correlation between the received patterns at different frequencies, which allows spectral information to be retrieved from multiplexed measurements. Frequency diversity is achieved by employing an aperture populated with metamaterial elements possessing diverse resonances. The frequency resolution in this implementation is determined approximately by the quality factor of these resonators [44]. Frequency correlations were computed as the equal-position complex product between different frequency components of the modulation functions, averaged over all positions and tuning states. That is, the frequency cross correlation is computed as

Correlation between our measurement modes is better quantified through the singular value decomposition, which gives an approximate measure of our ability to accurately reconstruct the target space from inversion of our linear model [59]. Figure 9(c) depicts the normalized singular value spectra computed for the matrix ${\bf H}$ as we increase the number of modulation patterns $M$. As the ratio of maximum to minimum singular values gives the condition number for the matrix, the nonzero slope observed in the singular value spectra indicates diminishing gains in the efficiency of our measured projections as correlations manifest in the random patterns. While more measurement modes are indeed accessed through the use of a larger number of modulation patterns, this number is effectively limited by noise. The tendency for measurement modes with small singular values to amplify noise in the inversion process motivates the use of a regularized solution for our experimental results.

Given an estimate of the CSD, the remaining characteristics of an image obtained with our system are determined by standard diffraction properties. In far-field imaging applications, a two-dimensional image is commonly achieved through a two-dimensional Fourier transform of the acquired coherence function, according to the Van Cittert–Zernike theorem [7]. In the near-field of the receiving aperture, a degree of longitudinal resolution is achievable due to the extent of longitudinal projections accessed by the aperture, enabling three-dimensional imaging through computational inversion via Eq. (12) [60,61]. This is evidenced by the compressive reconstruction of a single OEWG source in the aperture near-field at a position of $(30\;\text{cm}, - 2.8\;\text{cm},0.5\;\text{cm})$, shown in Fig. 10(a). Neglecting effects of the decoding procedure, system resolution can be approximately evaluated by examining the adjoint solution, i.e., the solution given by backprojecting a collection of standard array cross correlation measurements. We see from this case that, for an incoherent source, our imaging system obeys an incoherent imaging equation,

## 6. DISCUSSION

From the above results, we see that the proposed method can achieve passive, spectral imaging of sparse source distributions using a single DMA to encode the partially coherent fields onto single-port measurements. Several improvements in the system hardware and computational approach can yield more robust imaging performance. As indicated in Section 5, additional information can be gained simply through the use of more tuning states in the measurement and source estimation procedure up to the limit imposed by the aperture modulation capabilities, though this will result in longer measurement times. A more efficient measurement and estimation strategy would involve a more informed design of the aperture modulation functions, as opposed to the random modulation strategy employed here. For example, instead of random patterns comprising the encoding function ${A_m}(\vec r_a^\prime ,\omega)$, aperture modulation functions matched to the coherent modes [62] of a particular source could be formed dynamically by a properly designed aperture to achieve target-defined measurement acceleration.

Additional improvements in sensing performance can be achieved through more advanced coherence retrieval strategies. To this end, one may exploit physically informed priors on the coherence function solution [15,63] to obtain an accurate estimate of the full coherence function and complete characterization of a partially coherent system [64,65]. Under proper measurement calibration and synchronization conditions, multiple DMA panels can be integrated into a large-aperture system to improve image fidelity and resolution [58]. As our reconstruction model scales with scene size, the computational demands of matrix and iterative techniques can be cumbersome for large systems [66]. In such scenarios, alternative methods such as those based on Fourier domain techniques can offer significant reduction in computation time and power [26,60,67]. This interesting thrust is beyond the scope of this work and is pursued in future works.

We have restricted our studies in this work to scalar fields due to the use of a DMA design consisting of predominantly linearly polarized metamaterial elements. Using a metasurface with diverse polarization characteristics, one can potentially interrogate the full polarimetric structure of RF sources through a straightforward extension of the methods formulated here [68]. We also envision that, through proper inversion strategies mentioned above, the full coherence tensor of a partially coherent field can be retrieved using such polarization-diverse structures [69].

We have specifically investigated near-field passive imaging performance in order to recover three-dimensional, spectral microwave images. The problem of imaging sources in the aperture far-field is a simplified case of this technique, and can facilitate spectrally selective direction of arrival estimation for radar and multiple-input, multiple-output (MIMO) communication applications [70,71].

Finally, though the processing approach outlined here can be applied in principle to true thermal sources, the weak signal strength from such targets can pose a significant engineering challenge. A true thermal imaging system will likely require dedicated DMA and receiver design that properly optimizes metasurface antenna efficieny and receiver SNR.

## 7. CONCLUSION

In this work, we have formulated and demonstrated coded aperture passive microwave imaging as a method for spectral reconstruction of sources of incoherent radiation. We showed that this can be achieved using a single DMA operating in a lensless configuration using only spectral magnitude measurements, or total intensity measurements paired with suitable metasurface design. In contrast with alternative microwave array systems, our hardware is exceedingly simple and inexpensive, and the computational principles can be readily extended to terahertz and optical computational imaging systems.

## Funding

Air Force Office of Scientific Research (FA9550-18-1-0187).

## Disclosures

The authors declare that there are no conflicts of interest related to this paper.

## REFERENCES

**1. **J. W. Goodman, *Introduction to Fourier
Optics* (Roberts &
Company, 2005).

**2. **T. Pittman, Y. Shih, D. Strekalov, and A. Sergienko, “Optical imaging by means of
two-photon quantum entanglement,” Phys. Rev.
A **52**, R3429
(1995). [CrossRef]

**3. **J. H. Shapiro, “Computational ghost
imaging,” Phys. Rev. A **78**, 061802 (2008). [CrossRef]

**4. **A. V. Diebold, M. F. Imani, T. Sleasman, and D. R. Smith, “Phaseless computational ghost
imaging at microwave frequencies using a dynamic metasurface
aperture,” Appl. Opt. **57**, 2142–2149
(2018). [CrossRef]

**5. **A. V. Diebold, M. F. Imani, T. Sleasman, and D. R. Smith, “Phaseless coherent and
incoherent microwave ghost imaging with dynamic metasurface
apertures,” Optica **5**, 1529–1541
(2018). [CrossRef]

**6. **M. Soumekh, *Fourier Array Imaging*
(Prentice-Hall,
1994).

**7. **J. W. Goodman, *Statistical Optics*
(Wiley,
2015).

**8. **L. Mandel and E. Wolf, *Optical Coherence and Quantum
optics* (Cambridge
University, 1995).

**9. **K. Dutta and J. Goodman, “Reconstructions of images of
partially coherent objects from samples of mutual
intensity,” J. Opt. Soc. Am. **67**, 796–803
(1977). [CrossRef]

**10. **J. Garnier and G. Papanicolaou, *Passive Imaging with Ambient
Noise* (Cambridge University,
2016).

**11. **C. E. Yarman and B. Yazici, “Synthetic aperture hitchhiker
imaging,” IEEE Trans. Image Process. **17**, 2156–2173
(2008). [CrossRef]

**12. **S. J. Norton and M. Linzer, “Backprojection reconstruction
of random source distributions,” J. Acoust.
Soc. Am. **81**,
977–985 (1987). [CrossRef]

**13. **D. L. Marks, R. A. Stack, and D. J. Brady, “Astigmatic coherence sensor
for digital imaging,” Opt. Lett. **25**, 1726–1728
(2000). [CrossRef]

**14. **K. Nugent, “Wave field determination using
three-dimensional intensity information,”
Phys. Rev. Lett. **68**,
2261–2264 (1992). [CrossRef]

**15. **L. Tian, J. Lee, S. B. Oh, and G. Barbastathis, “Experimental compressive phase
space tomography,” Opt. Express **20**, 8296–8308
(2012). [CrossRef]

**16. **J. R. Fienup and P. S. Idell, “Imaging correlography with
sparse arrays of detectors,” Opt.
Eng. **27**, 279778
(1988). [CrossRef]

**17. **G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, “Generalized optical memory
effect,” Optica **4**,
886–892 (2017). [CrossRef]

**18. **J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non-invasive imaging through
opaque scattering layers,” Nature **491**, 232–234
(2012). [CrossRef]

**19. **A. R. Thompson, J. M. Moran, and G. W. Swenson, *Interferometry and Synthesis in Radio
Astronomy* (Wiley,
1986).

**20. **S. Kharkovsky and R. Zoughi, “Microwave and millimeter wave
nondestructive testing and evaluation - overview and recent
advances,” IEEE Instrum. Meas. Mag. **10**(2),
26–38 (2007). [CrossRef]

**21. **J. M. Beada’a, A. M. Abbosh, S. Mustafa, and D. Ireland, “Microwave system for head
imaging,” IEEE Trans. Instrum. Meas. **63**, 117–123
(2014). [CrossRef]

**22. **D. M. Sheen, D. L. McMakin, and T. E. Hall, “Three-dimensional
millimeter-wave imaging for concealed weapon
detection,” IEEE Trans. Microwave Theory
Tech. **49**,
1581–1592 (2001). [CrossRef]

**23. **M. Peichl, S. Dill, M. Jirousek, and H. Suess, “Passive microwave remote
sensing for security applications,” in
*European Radar Conference*
(IEEE, 2007),
pp. 32–35.

**24. **T. J. Jackson and T. J. Schmugge, “Passive microwave remote
sensing system for soil moisture: some supporting
research,” IEEE Trans. Geosci. Remote
Sens. **27**,
225–235 (1989). [CrossRef]

**25. **S. Vakalis and J. A. Nanzer, “Microwave imaging using noise
signals,” IEEE Trans. Microwave Theory
Tech. **66**,
5842–5851 (2018). [CrossRef]

**26. **J. Chen, Y. Li, J. Wang, Y. Li, and Y. Zhang, “An accurate imaging algorithm
for millimeter wave synthetic aperture imaging radiometer in
near-field,” Prog. Electromagn. Res. **141**, 517–535
(2013). [CrossRef]

**27. **C. S. Ruf, C. T. Swift, A. B. Tanner, and D. M. Le Vine, “Interferometric synthetic
aperture microwave radiometry for the remote sensing of the
earth,” IEEE Trans. Geosci. Remote
Sens. **26**,
597–611 (1988). [CrossRef]

**28. **N. Gopalsami, S. Liao, T. W. Elmer, E. R. Koehl, A. Heifetz, A. P. C. Raptis, L. Spinoulas, and A. Katsaggelos, “Passive millimeter-wave
imaging with compressive sensing,” Opt.
Eng. **51**, 091614
(2012). [CrossRef]

**29. **D. J. Brady, *Optical Imaging and
Spectroscopy* (Wiley,
2009).

**30. **M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel imaging via
compressive sampling,” IEEE Signal Process.
Mag. **25**(2),
83–91 (2008). [CrossRef]

**31. **G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded aperture
spectral imaging: an introduction,” IEEE
Signal Process. Mag. **31**(1),
105–115 (2013). [CrossRef]

**32. **J. Hunt, T. Driscoll, A. Mrozack, G. Lipworth, M. Reynolds, D. Brady, and D. R. Smith, “Metamaterial apertures for
computational imaging,” Science **339**, 310–313
(2013). [CrossRef]

**33. **M. F. Imani, J. N. Gollub, O. Yurduseven, A. V. Diebold, M. Boyarsky, T. Fromenteze, L. Pulido-Mancera, T. Sleasman, and D. R. Smith, “Review of metasurface antennas
for computational microwave imaging,” IEEE
Trans. Antennas Propag. **68**,
1860–1875 (2020). [CrossRef]

**34. **N. Landy, J. Hunt, and D. R. Smith, “Homogenization analysis of
complementary waveguide metamaterials,”
Photon. Nanostruct. Fundam. Applic. **11**, 453–467
(2013). [CrossRef]

**35. **T. Sleasman, M. F. Imani, J. N. Gollub, and D. R. Smith, “Dynamic metamaterial aperture
for microwave imaging,” Appl. Phys.
Lett. **107**, 204104
(2015). [CrossRef]

**36. **E. L. Kpré and C. Decroze, “Passive coding technique
applied to synthetic aperture interferometric
radiometer,” IEEE Geosci. Remote Sens.
Lett. **14**,
1193–1197 (2017). [CrossRef]

**37. **E. Kpre, C. Decroze, M. Mouhamadou, and T. Fromenteze, “Computational imaging for
compressive synthetic aperture interferometric
radiometer,” IEEE Trans. Antennas
Propag. **66**,
5546–5557 (2018). [CrossRef]

**38. **A. C. Tondo Yoya, B. Fuchs, and M. Davy, “Computational passive imaging
of thermal sources with a leaky chaotic cavity,”
Appl. Phys. Lett. **111**,
193501 (2017). [CrossRef]

**39. **A. Walther, “Radiometry and
coherence,” J. Opt. Soc. Am. **58**, 1256–1259
(1968). [CrossRef]

**40. **E. Wolf, “Coherence and
radiometry,” J. Opt. Soc. Am. **68**, 6–17
(1978). [CrossRef]

**41. **Y. Xie, T.-H. Tsai, A. Konneker, B.-I. Popa, D. J. Brady, and S. A. Cummer, “Single-sensor multispeaker
listening with acoustic metamaterials,” Proc.
Natl. Acad. Sci. USA **112**,
10595–10598 (2015). [CrossRef]

**42. **X. Li, J. A. Greenberg, and M. E. Gehm, “Single-shot multispectral
imaging through a thin scatterer,”
Optica **6**,
864–871 (2019). [CrossRef]

**43. **M. F. Imani, T. Sleasman, and D. R. Smith, “Two-dimensional dynamic
metasurface apertures for computational microwave
imaging,” IEEE Antennas Wireless Propag.
Lett. **17**,
2299–2303 (2018). [CrossRef]

**44. **T. Sleasman, M. F. Imani, A. V. Diebold, M. Boyarsky, K. P. Trofatter, and D. R. Smith, “Implementation and
characterization of a two-dimensional printed circuit dynamic
metasurface aperture for computational microwave
imaging,” arXiv:1911.08952v1
(2019).

**45. **I. Yoo, M. F. Imani, T. Sleasman, and D. R. Smith, “Efficient complementary
metamaterial element for waveguide-fed metasurface
antennas,” Opt. Express **24**, 28686–28692
(2016). [CrossRef]

**46. **B. T. Draine and P. J. Flatau, “Discrete-dipole approximation
for scattering calculations,” J. Opt. Soc. Am.
A **11**,
1491–1499 (1994). [CrossRef]

**47. **L. Pulido-Mancera, P. T. Bowen, M. F. Imani, N. Kundtz, and D. Smith, “Polarizability extraction of
complementary metamaterial elements in waveguides for aperture
modeling,” Phys. Rev. B **96**, 235402 (2017). [CrossRef]

**48. **L. Pulido-Mancera, M. F. Imani, P. T. Bowen, N. Kundtz, and D. R. Smith, “Analytical modeling of a
two-dimensional waveguide-fed metasurface,”
arXiv:1807.11592 (2018).

**49. **G. Lipworth, A. Rose, O. Yurduseven, V. R. Gowda, M. F. Imani, H. Odabasi, P. Trofatter, J. Gollub, and D. R. Smith, “Comprehensive simulation
platform for a metamaterial imaging system,”
Appl. Opt. **54**,
9343–9353 (2015). [CrossRef]

**50. **A. Yaghjian, “An overview of near-field
antenna measurements,” IEEE Trans. Antennas
Propag. **34**,
30–45 (1986). [CrossRef]

**51. **W. Menke, *Geophysical Data Analysis: Discrete
Inverse Theory* (Academic,
2018).

**52. **M. Elad, *Sparse and Redundant Representations:
From Theory to Applications in Signal and Image Processing*
(Springer,
2010).

**53. **J. M. Bioucas-Dias and M. A. Figueiredo, “A new twist: two-step
iterative shrinkage/thresholding algorithms for image
restoration,” IEEE Trans. Image
Process. **16**,
2992–3004 (2007). [CrossRef]

**54. **T. Fromenteze, O. Yurduseven, M. F. Imani, J. Gollub, C. Decroze, D. Carsenat, and D. R. Smith, “Computational imaging using a
mode-mixing cavity at microwave frequencies,”
Appl. Phys. Lett. **106**,
194104 (2015). [CrossRef]

**55. **O. Yurduseven, J. N. Gollub, D. L. Marks, and D. R. Smith, “Frequency-diverse microwave
imaging using planar Mills-Cross cavity apertures,”
Opt. Express **24**,
8907–8925 (2016). [CrossRef]

**56. **T. Shirai, T. Setälä, and A. T. Friberg, “Temporal ghost imaging with
classical non-stationary pulsed light,” J.
Opt. Soc. Am. B **27**,
2549–2555 (2010). [CrossRef]

**57. **M. F. Imani and D. R. Smith, “Temporal microwave ghost
imaging using a reconfigurable disordered cavity,”
Appl. Phys. Lett. **116**,
054102 (2020). [CrossRef]

**58. **J. Gollub, O. Yurduseven, K. Trofatter, D. Arnitz, M. Imani, T. Sleasman, M. Boyarsky, A. Rose, A. Pedross-Engel, T. H. Odabasi, G. L. Zvolensky, D. Brady, D. L. Marks, M. S. Reynolds, and D. R. Smith, “Large metasurface aperture for
millimeter wave computational imaging at the
human-scale,” Sci. Rep. **7**, 42650 (2017). [CrossRef]

**59. **D. L. Marks, J. Gollub, and D. R. Smith, “Spatially resolving antenna
arrays using frequency diversity,” J. Opt.
Soc. Am. A **33**,
899–912 (2016). [CrossRef]

**60. **D. L. Marks, R. A. Stack, and D. J. Brady, “Three-dimensional coherence
imaging in the Fresnel domain,” Appl.
Opt. **38**,
1332–1342 (1999). [CrossRef]

**61. **N. A. Salmon, “3-D radiometric aperture
synthesis imaging,” IEEE Trans. Microwave
Theory Tech. **63**,
3579–3587 (2015). [CrossRef]

**62. **A. Starikov and E. Wolf, “Coherent-mode representation
of Gaussian Schell-model sources and of their radiation
fields,” J. Opt. Soc. Am. **72**, 923–928
(1982). [CrossRef]

**63. **C. Bao, G. Barbastathis, H. Ji, Z. Shen, and Z. Zhang, “Coherence retrieval using
trace regularization,” SIAM J. Imaging
Sci. **11**,
679–706 (2018). [CrossRef]

**64. **H. N. Chapman, “Phase-retrieval x-ray
microscopy by Wigner-distribution deconvolution,”
Ultramicroscopy **66**,
153–172 (1996). [CrossRef]

**65. **J. N. Clark, X. Huang, R. J. Harder, and I. K. Robinson, “Dynamic imaging using
ptychography,” Phys. Rev. Lett. **112**, 113901
(2014). [CrossRef]

**66. **A. V. Diebold, L. Pulido-Mancera, T. Sleasman, M. Boyarsky, M. F. Imani, and D. R. Smith, “Generalized range migration
algorithm for synthetic aperture radar image reconstruction of
metasurface antenna measurements,” J. Opt.
Soc. Am. B **34**,
2610–2623 (2017). [CrossRef]

**67. **W. Carter, “On refocusing a radio
telescope to image sources in the near field of the antenna
array,” IEEE Trans. Antennas Propag. **37**, 314–319
(1989). [CrossRef]

**68. **T. Fromenteze, O. Yurduseven, M. Boyarsky, J. Gollub, D. L. Marks, and D. R. Smith, “Computational polarimetric
microwave imaging,” Opt. Express **25**, 27488–27505
(2017). [CrossRef]

**69. **J. Tervo, T. Setälä, and A. T. Friberg, “Theory of partially coherent
electromagnetic fields in the space–frequency domain,”
J. Opt. Soc. Am. A **21**,
2205–2215 (2004). [CrossRef]

**70. **D. D. Ross, C. J. Ryan, G. J. Schneider, J. Murakowski, and D. W. Prather, “Passive three-dimensional
spatial-spectral analysis based on k-space
tomography,” IEEE Photon. Technol.
Lett. **30**,
817–820 (2018). [CrossRef]

**71. **O. Yurduseven, M. A. B. Abbasi, T. Fromenteze, and V. Fusco, “Frequency-diverse
computational direction of arrival estimation
technique,” Sci. Rep. **9**, 1–12
(2019). [CrossRef]