Frequency Domain Image Reconstruction for Imaging With Multistatic Dynamic Metasurface Antennas

Dynamic metasurface antennas (DMAs) have recently been introduced as a computational imaging (CI) platform that offers significant advantages over conventional imaging systems. Such antennas are able to produce tailored radiation patterns that are later used to encode the scene’s information into a few measurements. However, since the signal is compressed, image reconstruction in the frequency domain using the conventional range migration algorithm (RMA) cannot be directly applied to CI-based apertures synthesized with DMAs. More sophisticated algorithms need to be developed to overcome this limitation, which transfers the complexity of such systems to the software layer. This paper proposes a new multistatic DMA RMA (MD-RMA) technique that achieves decompression of the compressed signals collected from multistatic DMA apertures and processes the decompressed data in the Fourier domain followed by a multistatic-to-monostatic conversion. Simulation results verify that the proposed approach can produce high quality 3D radar images, which alleviates the need to mechanically move the aperture. Moreover, since the image reconstruction is fully conducted in the frequency domain, leveraging the fast Fourier transform (FFT) algorithm, the complexity of the algorithm, and thus the execution time, is significantly reduced in contrast with conventional spatial domain algorithms.


I. INTRODUCTION
Microwave and millimeter wave imaging has been remarkably popular among a wide range of applications, such as security screening [1], through-wall imaging [2], and biomedical imaging [3]. Particularly in the field of security screening, the non-ionizing electromagnetic waves at those frequencies offer significant advantages [1], since they can penetrate common clothing without being hazardous for the human body. Although this fact brings microwave and millimeter imaging into the first line when developing such applications, there are The associate editor coordinating the review of this manuscript and approving it for publication was Mahmoud A. Abdalla . some challenges that security screening systems must meet to be successful. First, they need to provide the reconstructed images in real-time and, second, be cost-effective.
The traditional millimeter-wave imaging configurations are restricted to variants of synthetic aperture radar (SAR) [1], [4], [5], [6] or phased arrays systems [7], [8]. In the case of SAR, a co-located transmitting and receiving antenna is mechanically moved along one or two directions to synthesize an aperture. This raster scanning is complex in terms of hardware and is associated with slow acquisition time, especially when electrically large apertures are required, due to the mechanical movement of the antennas.
Phased arrays have eliminated the need of mechanical movement by establishing a stationary all-electronic system. However, this advantage does not counterbalance the many limitations of phased arrays. Such an aperture consists of an array, densely populated with antennas, which image the scene simultaneously. Considering the large apertures required for security screening applications, the number of antennas is large, leading to a complex system in terms of hardware. Moreover, each antenna needs extra electronic circuits to control its amplitude and phase, such as phase shifters and amplifiers, adding more complexity to the hardware. Considering the above, phased arrays have proven to be hardware complex, power hungry, and, as a result, an expensive alternative. These drawbacks have prevented the use of this technology in many applications.
Multiple-input multiple-output (MIMO) apertures have been also proposed in the literature [9], [10], [11], [12] to mitigate some of the above limitations. MIMO arrays can achieve all-electronic scanning, avoiding any mechanical movement, and employ fewer channels than phased arrays. Some hybrid configurations, called MIMO-SAR [13], [14], have also been developed to decrease mechanical scanning, while keeping the aperture size relatively small. Nevertheless, even in sparse scenarios [9], [12], the number of channels involved is notable.
In the last few years, computational imaging (CI) [15], [16], [17], [18], [19], [20], [21] has been proposed as an alternative to conventional systems. CI-based systems are based on antennas capable of producing quasi-orthogonal measurement modes that are able to acquire information of the investigation domain and compress it into a single (or a reduced number of) channels, rendering the data acquisition significantly faster. The probing of the investigation domain using the quasi-orthogonal measurement modes essentially eliminates the need for a raster-scan. In particular, metasurface antennas have proven to be more cost-effective, they are associated with reduced power consumption and are easy to fabricate [22]. On top of that, as required by CI, they can generate spatially incoherent, quasi-random radiation patterns to illuminate the whole scene at one time, in contrast to traditional SAR systems, for which mechanical movement is required. Also, compared to phased array systems, the overall complexity of the radiofrequency hardware is reduced as the additional circuitry needed to achieve the required phase-shifts is avoided.
The different measurement modes that are generated to encode the scene's information, namely, the quasi-orthogonal spatially incoherent radiation patterns, can be realized in different ways. One of them, called frequency diversity [18], [19], is based on the generation of different radiation patterns as a function of frequency. However, this requires the use of large bandwidths to acquire enough information. Another approach, which was recently introduced, is based on dynamic metasurface antennas (DMAs) [23], [24], [25], [26], [27], [28], [29]. A typical DMA consists of a waveguide with subwavelength metamaterial elements etched across its front surface. Utilizing an external voltage control, the elements can be switched on and off, such that a part of the elements radiates at each control level. Consequently, different on/off schemes are generated each time the voltage control changes, resulting in different radiation patterns that illuminate the scene. Each metasurface aperture pattern with a corresponding set of on/off metamaterial elements is called a mask. A sequence of radiation patterns is used to multiplex the scene information and encode it into a single or a reduced number of channels. As a result, compression of the physical layer is achieved.
Although a metasurface aperture mitigates the hardware constraints of conventional systems at a very satisfactory level, the price comes at the software level. As mentioned before, a DMA compresses the signal into a reduced number of channels, and, therefore, there is no longer a one-to-one relationship between each antenna element and each scene location. The scene information is associated with the compressed signal via a transfer matrix, usually called sensing matrix, and thus, spatial-domain CI algorithms mainly focus on solving an ill-conditioned inverse problem, using techniques such as matched filtering (MF) or least-squares (LS) [18], [19], [20], [21], [22], [24], [26], [30]. However, as the number of unknown variables increases, the scaling of these techniques becomes a bottleneck in terms of computational resources and time. As an example, image reconstruction using CI-based apertures and the conventional MF algorithm has been shown to require the processing of sensing matrices larger than 90 GB [19]. Clearly, this poses a significant burden on the signal processing layer of the imaging system.
Considering that for a wide range of applications, including security screening, the image reconstruction time is important, frequency domain image processing techniques [1] are preferable over spatial-domain CI methods, such as MF or LS [18], [19], [20], [21]. Frequency domain image reconstruction is based on Fourier transformations and can be realized using the fast Fourier transform (FFT) algorithm, decreasing significantly the computational complexity and thus, the reconstruction time. However, to apply the FFT algorithm, a uniform set of measurements is required, which is not the case when imaging with DMAs, since the signal is compressed. Hence, Fourier-based algorithms, like the range migration algorithm (RMA), cannot be directly applied to CI-based apertures synthesized with DMAs.
In [27], [28], [29], and [31], the authors demonstrate that if the field on the aperture is known, then the signal can be decompressed as if it had been acquired from a set of independent uniformly distributed magnetic dipoles along the DMA and propose a pre-processing step towards this direction. Then, backpropagation techniques, such as RMA, can normally be applied to obtain the reconstructed image. Nevertheless, the physical topology presented in those works is restricted to a multiple-input single-output architecture, consisting of one transmitting DMA along the horizontal axis and a single receiver probe to collect the backscattered infor-VOLUME 10, 2022 mation from the scene. This scenario is sufficient for twodimensional (2D) imaging across a single slice in the range plane. To obtain three-dimensional (3D) images, a mechanical movement of the DMA and/or the receiving probes is suggested. In [32], a panel-to-panel approach has been proposed, where one DMA serves as the transmitting aperture and another as the receiving one. In this implementation, a SAR scheme is adopted, where both DMA apertures are translated along one direction.
Motivated by increasing the receiving paths for the back-scattered signal and eliminating the need of timeconsuming mechanical scanning in case of 3D imaging, an all-electronic stationary multistatic structure based on DMA apertures is studied in this paper. From an image reconstruction perspective, a frequency domain processing technique, utilizing RMA applied to a multistatic DMA aperture (MD-RMA), is developed towards real-time operation of such a system. For this purpose, the idea of decompressing the signal presented in [27], [28], [29], and [31], is extended to be applicable to MIMO DMA configurations. The developed MD-RMA approach enables, for the first time, an FFT-based image reconstruction from CI MIMO apertures synthesized entirely by DMAs. A comparison between the proposed imaging technique and the techniques discussed above and provided in the references is presented in Table 1.
The rest of the paper is organized as follows: Section II introduces the mathematical model of a DMA and discusses the development of the pre-processing step required prior to applying the image reconstruction algorithm. The development of RMA is also addressed in this section. The ability of the proposed MIMO DMA-based configuration to provide high-quality 3D images is demonstrated through simulation results in Section III. Finally, Section IV summarizes the conclusions.

II. IMAGING WITH DYNAMIC METASURFACE ANTENNAS A. IMAGING CONFIGURATION
In this section, a multistatic scenario using DMAs for the transmitting and receiving apertures is demonstrated. The imaging configuration is depicted in Fig. 1. The transmitting DMAs (Tx), represented by the blue apertures, lie parallel to the x−axis, whereas the receiving DMAs (Rx), represented by the pink apertures, lie parallel to the y−axis.
In the case of Fig. 1, all apertures consist of the same number of metamaterial elements, which is denoted by N T for the transmitting DMA and N R for the receiving one, such that N T = N R . The center of the multistatic configuration is considered the center of coordinates. The center of the imaging scene lies z 0 meters away along z direction. For the depiction presented in Fig. 1, a single voxel scatterer lies in the center of the imaging scene.
The different tuning schemes, or masks, in the transmitting DMA are generated by changing the tuning voltage of the different metamaterial elements (they are switched on and off ). Each of these masks leads to a different radiation pattern in the transmitting aperture. Similarly, the receiving DMA cycles through different masks to receive the backscattered signal using different combinations of elements each time. The number of masks for the transmitting and receiving apertures is denoted by M T and M R , respectively, throughout this paper.
When using an antenna-to-antenna configuration, such as in conventional MIMO apertures [9], [10], [12], along with the first Born approximation [33], [34], the field in the aperture can be mathematically modeled using the formula: where σ ( p ) is the reflectivity of the scene and x T , y T , x R , y R , are the coordinates of each element on the transmitting and receiving apertures, respectively, and f is the frequency. In this case, the antennas are modeled like point-dipoles and the radiated field can be expressed using the Green's function as follows: where p refers to the coordinates of each element on the aperture, p to the coordinates of each voxel in the scene, and k is the corresponding wavenumber for each frequency f . However, in the configuration shown in Fig. 1, neither the receiver nor the transmitter can be modeled as single dipoles and a more complicated formula is required to mathematically describe the radiated fields on the aperture [27]. The measured signal can be written as: where m T and m R refer to the transmission and reception aperture masks respectively, while t and r refer to each transmission and reception DMA respectively. The terms U t,M T , U r,M R are used to describe the electric fields produced by each transmitting and receiving DMA [27].
Observing (3) and comparing it with (1), it is evident that there is no longer a straightforward relationship between the transmitting and receiving elements and the measured signal. In contrast, the signal in (3) is compressed into one channel for each Tx-Rx mask combination.

B. IMAGE RECONSTRUCTION
When imaging with the configuration in Fig. 1, the measured signal will be compressed as indicated by (3). To be able to use image processing techniques in the frequency domain, the signal must be written in the form of a typical multistatic signal as in (1). For this purpose, a pre-processing step needs to be developed to decompress the signal g as if it had been acquired from a uniform set of dipoles placed alongside each DMA. The former idea has been discussed in [27], [28], and [29], however, it needs to be redeveloped in order to be applicable for multistatic systems such as the one discussed in this paper.
The total fields over the transmitting and receiving apertures can be computed as follows [27]: where a m T /R (f ) models the behavior of each switched on metamaterial element. If the element is switched off for the specific mask, then the value of a m T /R (f ) is set to zero. Z 0 denotes the vacuum impedance, H 0 represents the guided magnetic field and β represents the propagation constant of the waveguide. Assuming that the signal measured on the aperture side can be expanded in terms of the fields involving all masks, it can be written as: where S(x i , t, y l , r, f ) is the signal measured from each independent dipole along the Tx, Rx DMAs. If the generated aperture modes experience some degree of orthonormality, then we can write: where + symbolizes the pseudo-inverse of a matrix, m T and m T , and m R andm R , identify the masks of the transmitting and receiving DMAs, respectively, and δ is the Dirac delta function.
Assuming that (7), (8) hold, then (6) can be inverted by multiplying both sides, first by the pseudo-inverse of the t,m T to the left and then by the pseudo-inverse of the r,m R to the right. After this step, we will have: or in matrix form: where T symbolizes the transpose of a matrix. According to (10), and considering that the aperture modes are chosen such that the required orthonormality is reassured, VOLUME 10, 2022 the data collected can be transformed as if they had been acquired by uniform distributed effective dipoles along each DMA. Since each matrix is extracted for each individual DMA, the transformation takes place for each combination of DMA transmitting-receiving apertures. The signal can then be written in the following form: where R T and R R are calculated as: respectively. Here, (x T , y T , z T ) and (x R , y R , z R ) represent the positions of Tx and Rx elements on each transmitting and receiving DMA, respectively, while (x , y , z ) denotes the scene's coordinates.
Considering that now the signal is expressed in terms of every transmitter and receiver location, it can undergo a multistatic-to-monostatic (MTM) transformation to obtain S(x c , y c , f ) [9], [10], [12]: whereŝ ref and s ref are reference signals and (x c , y c ) are the coordinates of the virtual grid. More details for their calculation can be found in [9], [10], and [12]. After this data conversion, the signal will lie on a uniform virtual grid (x c , y c ) on the original aperture plane. This grid is formulated based on the assumption that a phase-center is formed between each transmitter-receiver element, which constitutes the sampling position for the certain pair of elements in the far-field. Therefore, the signal S(x c , y c , f ) resembles the signal acquired from a monostatic SAR-like configuration [1] in which the sampling positions lie on the same grid, that is, (x c , y c ). It should be noted that this assumption is strictly valid in the far-field of the synthesized aperture. For near-field applications, a more sophisticated design of the aperture should be considered so that the effect of divergence from the far-field condition becomes negligible. More details are given in [9], [10], and [12], where multistatic-to-monostatic conversion has been successfully applied for near-field imaging scenarios.
Taking into account the above, the dispersion relationship is consequently defined similar to a conventional SAR scenario [1]: where k x , k y and k z represent the spatial wavenumber vectors.

C. ALGORITHM IMPLEMENTATION
The image reconstruction procedure is visualized in the diagram of Fig. 2. In this section, we dive into the implementation details for each algorithm step. The compressed signal measured in the aperture by a certain combination of Tx and Rx DMAs is denoted by g(t, r, f ) M T ×M R . The signal is then multiplied by + t (f ) N T ×M T and T + r (f ) M R ×N R from the left and right sides, respectively, for all combinations of transmitting and receiving DMAs, as indicated in (10). This constitutes the pre-processing step according to which the signal is decompressed as if it had been acquired from independent antennas along each DMA, so it can be written in the form of a typical multistatic signal, that is, S(x T , y T , x R , y R , k). This signal can be further processed using range migration techniques.
Next, the signal will undergo a multistatic-to-monostatic conversion as described in (14) to obtain the signal S(x c , y c , f ). Then a zero-padding operation takes place in the cross-range dimensions before the 2D Fourier transform. Such padding is essential before Fourier transformations to reduce aliasing in the frequency domain. The 2D FFT results inS(k x , k y , k), which is the signal in the frequency domain in terms of spatial wavenumbers k x , k y . Then the signal is back-propagated to the center of range direction (i.e. z− axis). This is achieved by multiplying by the term e jk z z 0 , where k z is also defined as indicated in (15). Vectors k x , k y are defined in the range [ −π d , π d ], where d denotes the sampling rate in the virtual grid (x c , y c ). Thus, d = λ min /2, i.e. half the wavelength at the maximum frequency. After this step we obtainS(k x , k y , k).
Then, the signal needs to be re-mapped in the k domain, such that it lies onto uniform points and a 3D IFFT can take place. This re-mapping is known as Stolt interpolation. Here, we are using a piecewise cubic hermitian interpolation (PCHIP), to address this step. After the interpolation, we obtainS(k x , k y , k z ), where k z consists of uniform samples and thus,S is in a compatible form for 3D IFFT. This is the final step of the reconstruction process, after which we obtain the final imageσ (x , y , z ).

III. SIMULATION RESULTS AND DISCUSSION
This section includes simulation results involving the setup shown in Fig. 1 and the adopted RMA algorithm described in the previous section. For all the simulation results presented in this paper, each DMA aperture resembles the one depicted in Fig. 1. The frequency ranges from 17.5 GHz to 22 GHz, such that the bandwidth is B = 4.5 GHz, sampled at n f = 32 equally distributed frequency points. The element spacing along the DMA is subject to a wavelength at 22 GHz. Consequently, the grid resulting from the multistatic-to-monostatic conversion, is sampled at half-wavelength at 22 GHz, [9], [10], [12], satisfying the Nyquist criterion. In this paper, N T = N R = 49 elements are considered on each DMA, resulting in apertures of physical size approximately equal to D = 65.5 cm. The number of masks is considered equal  to M T = M R = 40. For each mask, half of the elements will radiate while the rest remain switched off. The selection of the elements that are active for each mask is random. The imaging scene is the volume between x ∈ [−0.33, 0.33] m, y ∈ [−0.33, 0.33] m, z ∈ [1.6, 2.0] m. The distance from the synthesized aperture to the center of the imaging scene is z 0 = 1.8 m along the z direction.
As a first scenario, the imaging of a single point scatterer target, as depicted in Fig. 1, is carried out using the MD-RMA algorithm developed. The result of this study is important to analyze the point spread function (PSF) and confirm the diffraction-limited resolution of the synthesized multistatic DMA aperture. The reconstructed PSF pattern is illustrated in Fig. 3. In Fig. 4 the PSF curves are shown for both the  cross-range (xy− plane) and the range (zy− plane). The theoretical limits of the cross-range and range resolution for the results shown in Fig. 4 are: respectively, where λ min is the wavelength at maximum frequency, and c stands for the speed of light in vacuum. Given the simulation parameters described above, the resolution limits are computed, δ cr = 19 mm and δ r = 33 mm. Using the curves in Fig. 4, the cross-range resolution is calculated δ cr = 20 mm and the range one is δ r = 32 mm, confirming the above expectations. Following the investigation of the PSF characteristics of the synthesized multistatic DMA aperture, we consider the imaging of a distributed target consisting of a set of multiple point scatterers, as depicted in The reconstructed 3D image is shown in Fig. 6 as 2D cross sections across the xy−plane, xz−plane and yz−plane, respectively. Specifically, Figs. 6(a), 6(d), and 6(g) illustrate the reconstructed image when the proposed MD-RMA technique is used. The dynamic range for all figures is equal to 5 dB, and the side colorbar refers to the distance (in meters) of each target from the aperture. As can be seen in Fig. 6, the point scatterers are clearly retrieved in the reconstructed images in both the cross-range and range planes. Fig. 6 also provides a comparison between the proposed MD-RMA reconstruction algorithm and the MF and LS algorithms, which address the problem in spatial domain [18], [19], [20], [21]. The MF technique estimates the reflectivity of the scene as follows: where H represents the sensing matrix and † the conjugate transpose [18], [19], [20], [21]. Figs. 6(b), 6(e), and 6(h) present the 3D reconstructed image of the same target when MF is used. On the other hand, LS is an iterative technique that tries to converge to a better solution through a non-linear process. In particular, it aims to minimize the following cost function: The solution column vector derived form (18) was used as initial guess for the iterative method. Figs. 6(c), 6(f), and 6(i) present the 3D reconstructed image of the same target when the LS method is used. The scene is discretized into 17 × 17 × 31 voxels, i.e. 8959 voxels in total. In Figs. 6(b), 6(e), and 6(h), we observe that in the case of MF, the resolution is better for targets located closer to the aperture. However, in the MD-RMA results ( Fig. 6(a,d,g)), the resolution improves for the targets placed further from the aperture. This is because, in the proposed algorithm, the phase-center approximation has been considered. This approximation is more accurate as we move towards the far-field of the aperture, especially for corner-point scatterers away from the optical axis [9], [12]. Finally, in the case of LS method (Fig. 6(c,f,i)), the resolution is improved compared to the other two methods and regardless of the distance from the targets to the aperture. However, this method requires significant computational resources and might have convergence problems in the case of noisy measurements.
Simulation results were obtained using MATLAB 2021b on a computer with 64-bit Windows 10 operating system, 16GB of random access memory (RAM), and a Core i7-8700 central processing unit (CPU). The execution time for the MD-RMA algorithm of Fig. 2 is included and compared to the MF and LS techniques in Table 2. The times shown are the average time of ten consecutive executions on the CPU. As decompression, we consider the time needed to express the signal in terms of independent dipoles along each DMA, i.e. executing (10). By MTM conversion, we refer to the part of the algorithm that transforms the multistatic signal to lie on the virtual uniform grid, that is, the step described by (14). By RMA we refer to the frequency domain part of the algorithm, i.e. from applying the 2D FFT to the final 3D image, as seen in flow chart in Fig. 2. The total execution time of the proposed approach is the sum of the aforementioned parts and is shown in the last row of the Table 2. For the MF and LS approaches, we only need to consider the time needed to perform (18) and (19) respectively, considering that the sensing matrix is pre-computed since it is not target dependent. The last row of the second and third columns refers to this execution time. Looking at the bottom row of Table 2, a significant advantage of the proposed MD-RMA algorithm over the MF approach can be appreciated in terms of the execution time. It is evident that the proposed MD-RMA approach is executed much faster on the CPU, in the order of 30x, including all pre-processing operations needed. Taking into account that the LS algorithm utilizes the reconstruction of MF as initial guess and solves (19) in an iterative manner, its execution time is even larger. From Table 2, it can be concluded that the execution time of the proposed MD-RMA algorithm is only 0.32% of the time needed to execute the LS algorithm.

IV. CONCLUSION
In this paper, a stationary multistatic DMA-based CI system was proposed for 3D imaging at microwave frequencies. The proposed system is capable of producing good quality images, alleviating the need for mechanical movement and the use of numerous channels. To achieve image reconstruction for the synthesized CI-based DMA aperture, we developed a multistatic RMA technique, called MD-RMA, compatible with the presented multistatic DMA-based system. The execution time of the overall reconstruction procedure is promising (on the order of seconds) and exhibits a substantial improvement in comparison to the MF and LS algorithms. Employing general purpose graphical processing units (GPGPUs), this time can be further reduced, making a step towards real-time imaging with DMA-based MIMO systems. Her current research interests include inverse scattering, remote sensing, radar systems, imaging techniques, antenna measurement and diagnostics, and non-invasive measurement systems on board unmanned aerial vehicles. She has authored more than 40 peer-reviewed journals and conference papers, and she holds two patents. She has also received several awards, such as the 2020 National Award to the Best Ph.D. Thesis on telecommunication engineering.
GUILLERMO ÁLVAREZ-NARCIANDI received the M.Sc. degree in telecommunication engineering and the Ph.D. degree from the University of Oviedo, Gijón, Spain, in 2016 and 2020, respectively. He was a Visiting Student with Stanford University, Stanford, CA, USA, in 2014, and a Visiting Scholar with the University of Pisa, Pisa, Italy, in 2018, and the Institute of Electronics, Microelectronics and Nanotechnology, University of Lille, Lille, France, in 2019. In 2022, he joined Queen's University Belfast, Belfast, U.K., as a Research Fellow with the Centre for Wireless Innovation. His research interests include radar systems and imaging techniques, antenna diagnosis and characterization systems, localization and attitude estimation systems, and RFID technology. He has received the AMTA 2019 Student Paper Award (Second Place) and the Special Award to the Best Entrepreneurship Initiative in XV Arquímedes National Contest, in 2017, for the development of a RFID-based location system.
OKAN YURDUSEVEN (Senior Member, IEEE) received the Ph.D. degree in electrical engineering from Northumbria University, Newcastle upon Tyne, U.K., in 2014. He is currently a Reader (an Associate Professor) at the School of Electronics, Electrical Engineering and Computer Science, Queen's University Belfast, U.K. Prior to this, he was a NASA Research Fellow at the NASA Jet Propulsion Laboratory, California Institute of Technology, USA. He is also an Adjunct Professor at Duke University, USA. His research interests include microwave and millimeter-wave imaging, multiple-input-multiple-output (MIMO) radars, wireless power transfer, antennas and propagation, and metamaterials. He has authored more than 150 peer-reviewed technical journals and conference papers. He has been the principal investigator and the co-investigator on research grants totaling in excess of £10M in these fields.