Trading beams for bandwidth: Imaging with randomized beamforming

TL;DR

This paper introduces a broadband imaging technique that reduces the number of spatial measurements needed by using randomized beamforming, leveraging the scene's range-limited structure to improve efficiency and cost-effectiveness.

Contribution

It proposes a novel trade-off allowing fewer array measurements through broadband illumination and randomized beamforming, supported by theoretical analysis and simulations.

Findings

Significant reduction in spatial measurements needed for imaging

Theoretical validation of the measurement trade-off

Improved imaging speed and cost efficiency

Abstract

We study the problem of actively imaging a range-limited far-field scene using an antenna array. We describe how the range limit imposes structure in the measurements across multiple wavelengths. This structure allows us to introduce a novel trade-off: the number of spatial array measurements (i.e., beams that have to be formed) can be reduced significantly lower than the number array elements if the scene is illuminated with a broadband source. To take advantage of this trade-off, we use a small number of "generic" linear combinations of the array outputs, instead of the phase offsets used in conventional beamforming. We provide theoretical justification for the proposed trade-off without making any strong structural assumptions on the target scene (such as sparsity) except that it is range limited. In proving our theoretical results, we take inspiration from the sketching literature. We also provide simulation results to establish the merit of the proposed signal acquisition strategy. Our proposed method results in a reduction in the number of required spatial measurements in an array imaging system and hence can directly impact their speed and cost of operation.

Tables

Table 1. Table 1: Relative reconstruction error values for different classes of images. Aperture codes with weights chosen from standard normal distribution were used. The last column shows the output SNR as calculated using equation ( 31 ).

Imaging mode	Constant range (CR)	Flat surface (FS)	Multi-depth, known range profile (MD)	Multi-depth, unknown range profile (MDUR)	Constant range noisy
					Signal error	OP SNR
Full (1100 beams)	NA	NA	NA	NA	NA	16.09
320 codes	2.7e-5	4.4e-4	3e-4	1.1e-2	2.5e-5	16.0117
160 codes	7.4e-5	3.3e-4	1.2e-3	3.1e-2	6.3e-5	15.9112
80 codes	4.2e-4	5.1e-4	4.8e-2	1.6e-1	4.9e-4	15.3206

Table 2. Table 2: Relative reconstruction errors with subsampled array: The array was randomly subsampled to have 320, 160 and 80 elements. Abbreviations are as in Table 8 .

Imaging mode	CR	FS	MD	MDUR
320 elements	2e-4	1e-4	1e-3	3.8e-2
160 elements	9e-4	4e-4	3.1e-2	6.8e-2
80 elements	1.6e-3	1.1e-3	5.8e-3	18.4e-2

Equations

$$y_{\lambda }=A_{\lambda }{x}_0$$

$$z_{\lambda }={\psi }_{\lambda }{y}_{\lambda }={\psi }_{\lambda }{A}_{\lambda }{x}_0$$

$$\underset{x}{\min }\Vert \Phi y-\Phi Ax{\Vert }^2$$

$$\underset{x}{\min }\Vert y-Ax{\Vert }^2$$

$${\tilde{y}}_{d,\lambda }(t)={\mathrm{e}}^{\mathrm{j}2\pi ct/\lambda }\underset{-\pi /2}{\overset{\pi /2}{\int }}\int p(r,\theta ){\mathrm{e}}^{-\mathrm{j}2\pi (r+d\sin \theta )/\lambda }drd\theta .$$

$$y_{d,\lambda }=\underset{-\pi /2}{\overset{\pi /2}{\int }}\int p(r,\theta ){\mathrm{e}}^{-\mathrm{j}2\pi (r+d\sin \theta )/\lambda }drd\theta ={\hat{x}}_c\left({\omega }_r,{\omega }_{\tau }\right).$$

$${\tilde{y}}_d(t)=\underset{-\pi /2}{\overset{\pi /2}{\int }}\int p(r,\theta )s_b(t-r/c-d\sin \theta )drd\theta .$$

$$y_c=\mathcal{S}\mathcal{F}{x}_c=\{{\hat{x}}_c({\omega }_{\tau },{\omega }_r):({\omega }_{\tau },{\omega }_r)\in \mathcal{W}\}$$

$${\tilde{z}}_{i,\lambda }(t)=\sum _{m=1}^M{\varphi }_{\lambda }(i,m){\hat{y}}_{m,\lambda }(t)$$

$$z_{\lambda }={\varphi }_{\lambda }{y}_{\lambda }$$

$$z_{\lambda }=\varphi y_{\lambda }$$

$$y_c=\mathcal{S}\mathcal{FQ}x.$$

$$y_m(\lambda )$$

$${\tilde{y}}_{m,\lambda }(t)={\mathrm{e}}^{\mathrm{j}2\pi (ct+R_0)/\lambda }\underset{-1/2}{\overset{1/2}{\int }}p(R_0,\tau ){\mathrm{e}}^{-\mathrm{j}2\pi {\omega }_{\tau }\tau }\mathrm{d}\tau .$$

$$y_m(\lambda )=\underset{-1/2}{\overset{1/2}{\int }}p(R_0,\tau ){\mathrm{e}}^{-\mathrm{j}2\pi 2d_m\tau /\lambda }\mathrm{d}\tau .$$

$$g({\omega }_{\tau })=\underset{-1/2}{\overset{1/2}{\int }}p(R_0,\tau ){\mathrm{e}}^{-\mathrm{j}2\pi {\omega }_{\tau }\tau }\mathrm{d}\tau .$$

$$A_{R_0,\lambda }(m,n)={\exp }^{-j2\pi R_0/\lambda }{\exp }^{-i2\pi 2d_m(-0.5+n/N)/\lambda }.$$

$$y_{R_0,\lambda }=A_{R_0,\lambda }{x}_{R_0}.$$

$$A_{R,\lambda }(m,n)={\exp }^{-j2\pi (2d_m(-0.5+n_{\tau }/N)+(R_{\min }+d_{r}R/D))/\lambda }$$

$$y_{R,\lambda }=A_{R,\lambda }{x}_R.$$

$$y=\left[\begin{array}{c}{y}_{{\lambda }_1}\\ y_{{\lambda }_2}\\ ⋮\\ y_{{\lambda }_K}\end{array}\right]=\left[\begin{array}{c}{A}_{{\lambda }_1}\\ A_{{\lambda }_2}\\ ⋮\\ A_{{\lambda }_K}\end{array}\right]x_0=A{x}_0.$$

$$z=\Phi y=\left[\begin{array}{cccc}\varphi & 0& \cdots & 0\\ 0& \varphi & \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \varphi \end{array}\right]\left[\begin{array}{c}{y}_{{\lambda }_1}\\ y_{{\lambda }_2}\\ ⋮\\ y_{{\lambda }_K}\end{array}\right]=\Phi A{x}_0.$$

$$\underset{x}{arg\; min}\Vert \Phi y-\Phi Ax{\Vert }_2^2.$$

$$\underset{x}{arg\; min}\Vert y-Ax{\Vert }_2^2.$$

$$x_{\text{SLS}}=\text{arg}\underset{x}{\min }\Vert Sy-SAx{\Vert }^2,$$

$$\Phi y=\Phi A{x}_0=Y{x}_0,$$

$$\Phi =\left[\begin{array}{cccc}\varphi & 0& \cdots & 0\\ 0& \varphi & \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \varphi \end{array}\right].$$

$$\left[\begin{array}{c}{A}_1\\ A_2\\ \cdot \\ A_k\end{array}\right]=\left[\begin{array}{cccc}{C}_{11}& 0& \cdots & 0\\ C_{21}& C_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C_{k1}& C_{k2}& \cdots & C_{kk}\end{array}\right]\left[\begin{array}{c}{V}_1^T\\ V_2^T\\ ⋮\\ V_k^T\end{array}\right]=C{V}^T$$

$$x_{\text{LS}}=A^{†}y=A^{†}A{x}_0=V{V}^{\ast }{x}_0$$

$$V_{\Phi A}{V}_{\Phi A}^{\ast }-V{V}^{\ast }=0.$$

Full text

1 Introduction

In active array imaging, a transmitter emits an excitation signal and then forms an image using the reflections collected at an array of sensors. This technique has been employed in a multitude of fields ranging from medicine to security and surveillance, among many others. Array imaging offers a window outside of the visible spectrum, which can prove crucial in many applications where visible light cannot penetrate. However, a high cost barrier has historically prevented widespread adoption of array imaging in commercial products [1]. Nevertheless, applications such as autonomous vehicles, depth sensing, gesture recognition [2] and others have caused an increase of interest in commercializing active imaging modalities such as LiDAR and RADAR, leading to a number of efforts aimed at reducing the cost and increasing the efficiency of array imaging systems. In this paper, we consider a general antenna array imaging system and propose a novel trade-off that utilizes bandwidth of the excitation signal to reduce the number of array measurements needed to image targets that are range-limited, thus directly impacting the speed and cost of acquisition.

††This work was supported, in part, by NSF grants CCF-1350616 and CCF-1409406, a grant from Lockheed Martin, and a gift from the Alfred P. Sloan Foundation††A part of the work presented here appeared in [3] and in [4]

We use the standard linear model for array imaging in a setting with a single transmitter and multiple receiver elements. We use $M$ to denote the number of antenna array elements and $\lambda$ to denote the wavelength of the active excitation signal. We will later show that the set of discrete measurements collected by the antenna array elements at wavelength $\lambda$, $y_{\lambda}$, can be modeled as

[TABLE]

where $x_{0}\in\mathbb{R}^{N}$ denotes the sampled target scene and $A_{\lambda}$ is the linear operator ($M\times N$ matrix) mapping the scene to the measurements. Collecting broadband measurements is equivalent to collecting $y=\{y_{\lambda}\}$ for ${\lambda=\lambda_{1},\ \lambda_{2},\ \cdots,\ \lambda_{K}}$. We note that all $\{y_{\lambda}\}$ can be obtained simultaneously (using only $M$ “spatial measurements”) by using a single broadband excitation signal and computing a temporal Fourier transform.

A standard imaging method known as beamforming collects linear combinations of the array elements’ outputs instead of sampling them directly. This can be modeled as:

[TABLE]

where $\psi_{\lambda}$ is typically an invertible $M\times M$ matrix. In traditional beamforming the weights are chosen to induce spatial selectivity where each of the $M$ measurements collect reflections from distinct spatial sectors. This can be considered a special case of aperture coding: collecting linear combinations of the array outputs. Henceforth, we refer to linear combinations of array outputs as spatial measurements or beams. Note that beams conventionally refer to directivity inducing linear combinations, but for the purposes of this paper, we refer to any linear combination of the array outputs as a beam.

Our main observation is that the set of spatial, broadband measurements of such range-limited targets lie in a lower dimensional subspace and hence have a limited number of degrees of freedom. This naturally leads to the question of whether spatially undersampled array data of such scenes can still yield reconstructions as good as those obtained with full data. We provide a positive answer to this question. While this is reminiscent of compressed sensing, an important distinction between our work and the compressed sensing paradigm is that we do not impose any model such as sparsity on the target scene. We only require the scene to be range-limited.

Our motivation in this work is based on the observation that sampling and storing all of the array elements’ outputs, or obtaining all of the $M$ possible linear combinations required for traditional beamforming can be challenging and is in fact wasteful when the target scene is range-limited. We show that in this case, by taking fewer generic linear combinations (or aperture codes) or even just spatial subsampling, one can obtain the same quality reconstructions as that of using full measurements.

We take inspiration from a set of dimensionality reduction techniques used in the numerical linear algebra community known as sketching [5] to provide theoretical justification for the proposed signal acquisition method. We show that image reconstruction with a few generic aperture codes is equivalent to a sketched least squares problem of the form

[TABLE]

where $\Phi$ is a highly structured compresssive matrix. We show that it provides a solution equivalent to that of solving the higher dimensional problem

[TABLE]

which corresponds to standard image reconstruction methods. However, we again emphasize that unlike compressed sensing techniques, we do not assume any sparsity in the image domain.

As an example, consider an imaging setup where the array is two dimensional and has $40\times 40$ sensors with sensors placed $\lambda/2=3.75$cm apart and a target scene having a span of $[-45^{\circ},45^{\circ}]$ in both elevation and azimuthal angles. Let the scene have delta thickness: only one reflector per each angle, present at a constant known depth. Standard ways of imaging such a scene would require around $1100$ beams at wavelength $\lambda=7.5$cm. By introducing bandwidth in the excitation signal, we show that the scene can be imaged with as few as $80$ spatial linear combinations. This is illustrated in Figure 1. For target scenes with higher range limits, standard imaging methods also need bandwidth for imaging [6]. We however show that this bandwidth can be used to obtain similar gains in the number of spatial linear combinations and provide theoretical justification for the gains that aperture coding/subsampling can provide.

2 Related work

Array imaging has been addressed in a number of previous works. Standard array imaging problems have been considered in [6, 7] where the inverse problem is set up using transmit and receive beamforming with narrowband excitation. The authors also describe the fundamental limits on the resolution of array imaging systems arising due to finite apertures and discrete sampling of the array. General 3D imaging with wide-band excitation is considered in [8]. In general, to identify a 3D scene, a 2D antenna array and wideband excitation are necessary. Our focus is on a different dimension of this problem: we identify the limited number of degrees of freedom in a range-limited image and explore an alternate way using fewer measurements in which such a scene can be captured.

As arrays get larger, or when arrays need to be used in low-cost commercial applications, it is desirable to reduce the reduce the number of elements in an array and to use smarter algorithms to reduce the cost of the system. By sequentially using different parts of the array, one can obtain a set of low resolution images and propose to achieve reconstruction by upsampling and summing these images [9, 10]. In another approach, using a carefully designed non-uniform array, [11] proposes to increase the number of resolvable directions to $O(N^{2})$ by using an array with $O(N)$ elements. This approach enables the use of specially designed arrays to solve the problem of direction of arrival of source signals in passive sensing scenarios. Another main theme in reducing the cost of array systems has been the use of compressed sensing techniques. Reducing sampling rates at the sensors for digital beamforming is proposed in [12, 13, 14], but the number of beams/ array elements remain the same as conventional imaging. In a slightly different application, [15] imposes sparsity based regularization to solve the ill-posed radio-interferometric imaging problem. A similar theme is also followed in [16] which addresses the same problem of 3D imaging with a 2D phased array, but assumes sparsity in the image domain. The trade-off we propose has a different flavor: we demonstrate that range-limited images have a limited number of degrees of freedom and impose no further models such as sparsity or total variation on the image. Unlike in [14, 13, 12], we directly address the number of measurements/array elements rather than sampling rates at each element.

The problem of understanding the degrees of freedom in the context of active imaging of range-limited scenes is a particular case of the phenomenon of simultaneous concentration of energy of a signal in spatial/temporal and spectral domains. Signals which are concentrated maximally in a given time interval and a frequency interval are well-approximated using a subspace of dimension approximately equal to the product of the lengths of the intervals. This has been well studied for one dimensional signals [17, 18, 19, 20, 21, 22, 23]. Spatio-spectral concentration in two dimensions has also been studied in [24, 25]. We study this phenomenon in the case antenna array imaging and propose methods to drive the number of measurements close to the number of degrees of freedom, in contrast to conventional imaging methods that collect far more number of measurements.

We model our approach as a novel matrix sketching problem and provide guarantees on the number of measurements/beams required to image a far-field target. Matrix sketching refers to a set of techniques in numerical linear algebra for dimensionality reduction. In particular, a given large matrix is pre- or post-multiplied by a suitable randomized matrix to reduce the ambient dimensions while still retaining the required information, such as the Euclidean distance between elements in its column space. Research in this area has seen a rise in popularity due to its utility in solving large problems in numerical linear algebra [26, 5]. Our theoretical results are most closely aligned with those in [27], where random Gaussian projections are used to capture the range space of linear operators. In the context of active imaging, usage of sketching ideas can be seen in [28] where the Fourier basis is used as a sketching matrix for dimensionality reduction in interferometry. But the sketching matrix in [28] is a generic rectangular sketching matrix. In the context of wideband array imaging, the sketching matrix involved itself has a very particular structure that is dictated by the physical setup, which introduces new challenges. The work presented in the remainder of this paper is an extension of ideas first presented in a preliminary form in [4, 3].

3 Active array imaging

In this section, we describe the standard setup and measurement models of active array imaging. While doing so, we highlight some facts about broadband array imaging that form the basis of our contribution. In particular, we will emphasize that broadband array measurements provide a set of bandlimited Fourier domain samples of the target scene. We also describe some standard acquisition methods used to physically collect these samples. In the subsequent sections, we show that these bandlimited samples are highly structured for range-limited target scenes.

3.1 Propagation model and Fourier domain samples

We first consider the broadband array imaging problem with a one-dimensional (1D) antenna array. This can be easily extended to two dimensions. Consider a uniform linear array of aperture length $D$ placed on the $x$-axis, spanning $[-D/2,D/2]$. For now, we assume that the aperture is continuous. The phase center is taken at the origin: i.e., the time delay of arrival at each element is measured with respect to the element at the origin. We assume that the scene to be imaged lies in the X-Y plane and is in the far-field region of the array. Imaging the scene is equivalent to reconstructing its reflectivity map, which is a function of the roundtrip distance of the target from the array center and the angle from the broadside, to be denoted as $p(r,\theta)$. The system consists of only one transmitting element which is co-located with the receiver at the array center. Let $t$ be the continuous time index. For an excitation signal $s(t)$, the signal received at the array location $d$ is $s(t-r_{0}/c-d\sin\theta_{0}/c)$ for a unit-strength reflector at $(r_{0},\theta_{0})$. For a general reflectivity map, the narrowband response at this location for the excitation signal $s(t)=e^{j2\pi ct/\lambda}$ is

[TABLE]

By making the substitution $\tau=(\sin\theta)/2$, the complex amplitude of the signal received at location $d\in[-D/2,D/2]$ for excitation wavelength $\lambda$ can be written as

[TABLE]

where $\omega_{r}=\frac{1}{\lambda},\ \omega_{\tau}=\frac{2d}{\lambda}$, and $\widehat{x}_{c}$ denotes the Fourier transform of $x_{c}(r,\tau)=\frac{p(r,\sin^{-1}(2\tau))}{\sqrt{1-4\tau^{2}}}$. This shows that the antenna aperture measures the Fourier transform of (a warped version of) the target scene. $y_{d,\lambda}$ can be obtained by computing the Fourier transform of the temporal signal $\widetilde{y}_{d,\lambda}(t)$ received at the array after sampling it at a suitable rate, or by measuring the complex amplitude of the received signal. As a direct consequence of the finiteness of the aperture,we have that at any excitation wavelength $\lambda$, the accessible interval in the Fourier domain for $\omega_{\tau}$ is limited to $[-D/\lambda,D/\lambda]$.

For an imaging system with a 2D array in the X-Y plane and a 3D scene, the extension of the setup is straightforward to derive. The coordinates in 3D can be denoted as $(r,\theta_{1},\theta_{2})$, where $r$ is the roundtrip distance to the array center, $\theta_{1}$ is the angle with respect to the Y-Z plane, and $\theta_{2}$ is the angle with respect to the X-Z plane, as shown in Figure 2b. The scene reflectivity is denoted as $x_{c}(r,\tau_{1},\tau_{2})$ where $\tau_{1}=(\sin\theta_{1})/2$ and $\tau_{2}=(\sin\theta_{2})/2$. At excitation wavelength $\lambda$, the 2D array outputs $y_{d_{1},d_{2},\lambda}=\widehat{x}_{c}(1/\lambda,\,2d_{1}/\lambda,\,2d_{2}/\lambda)$ are samples of the Fourier transform of the scene sampled in the region bounded by \big{(}\frac{\pm D_{1}}{\lambda},\frac{\pm D_{2}}{\lambda}\big{)}, where $D_{1}$ and $D_{2}$ are the dimensions of the 2D array. From now on, we use the 1D array to discuss our model for the sake of notational brevity. However, all our discussion and results hold for both cases and all simulations use 2D arrays.

We are mainly interested in a broadband excitation scenario. We assume that the excitation signal is a broadband pulse with bandwidth in the interval $[\lambda_{\min},\lambda_{\max}]$. If the broadband signal used is $s_{b}(t)$, then the received signal at location $d$ is

[TABLE]

The complex amplitudes at different wavelengths, $y_{d,\lambda}$ can then be obtained as the Fourier coefficients of $\widehat{y}_{d}$ after taking its temporal Fourier transform. Since the lateral frequency support increases with decrease in the excitation wavelength, the accessible Fourier domain has a trapezoidal shape as illustrated in Figure LABEL:sub@fig:wedge_a. Similarly, for a 2D antenna, the region in the Fourier domain measured is shown in Figure LABEL:sub@fig:wedge_b.

It is now clear that a finite aperture and a finite bandwidth of excitation together offer bandlimited measurements of the target scene. If $\mathcal{F}$ represents the continuous domain Fourier transform and $\mathcal{S}$ represents the bandlimiting operator supported only in this accessible window, the continuous aperture measurements $y_{c}=\{y_{d,\lambda}\}$ can be modeled as

[TABLE]

where $\mathcal{W}$ represents the trapezoidal region shown in Figure 3. Broadband imaging is hence the task of collecting Fourier measurements in a bandlimited region and inferring the target reflectivity profile using these measurements.

3.2 Array measurement model

Our discussion above describes what can be observed through a finite aperture. In practice, we must measure this signal using a discrete array of sensors. This limits the Fourier domain measurements considered in the previous section to only a discrete set of samples. Let the 1D antenna considered in the previous section be an array of $M$ discrete antenna elements placed uniformly at coordinates $d_{\frac{-M}{2}},\cdots,d_{\frac{M}{2}}\in[-D/2,D/2]$. For ease of explanation, we consider a discrete set of $K$ excitation wavelengths $\{\lambda_{1}=\lambda_{\min},\lambda_{2},\cdots,\lambda_{K}=\lambda_{\max}\}$. When the complex amplitudes are measured at these wavelengths, the Fourier samples obtained are located on a pseudopolar grid, as shown in Figure 4. Measurements at these wavelengths can be obtained by using a single broadband excitation, as explained earlier. The set of all measurements $\{y_{m,\lambda}\}$ where $m=1,\cdots,M$, $\lambda\in\{\lambda_{1},\cdots,\lambda_{K}\}$ can be denoted by a vector $y\in\mathbb{C}^{MK}$. The set of measurements at each wavelength $\lambda$ is denoted by $y_{\lambda}\in\mathbb{C}^{M}$.

Traditional imaging involves collecting the samples shown in Figure 4 to reconstruct the target reflectivity profile. This can be achieved in many ways. One way is to directly read out the output of each antenna element. This amounts to measuring $y$ directly. Physically, this can be realized by using a broadband pulse as before and then taking a Fourier transform. Alternatively, one could cycle through a set of narrowband excitation signals (stepped frequency excitation) and collect array measurements at each wavelength. In either case, this method would require $M$ array element readouts.

Another standard method for acquiring the measurements is to collect linear combinations of the array outputs. In this method, the output of each array element is weighted and then added to the output of other elements. In particular, the procedure known as beamforming obtains specific linear combinations of the array outputs that induce spatial directivity. In narrowband beamforming, the weights are chosen such that the time delays for signals coming from a particular spatial direction are compensated for, hence “focusing” the array in that physical direction, as outlined in [6]. Instead of acquiring $M$ direct read-outs, $M$ linear combinations are acquired.

Acquiring linear combinations at each excitation wavelength separately is equivalent to acquiring linear combinations of samples along each row of points in Figure 4. A different set of linear combinations may be used for each excitation wavelength. Mathematically, if the excitation wavelength is $\lambda$, the measurements made in time domain are

[TABLE]

for $i=1,\cdots,M$. The vector of complex amplitudes $z_{\lambda}\in\mathbb{C}^{M}$ is then given by

[TABLE]

where $\phi_{\lambda}$ is in general an $M\times M$ well-conditioned matrix whose $(i,m)^{\text{th}}$ element is $\phi_{\lambda}(i,m)$. If a single broadband pulse is used for excitation, then the set of weights for the linear combinations at different wavelengths are constrained to be the same and the vector of complex amplitudes at different wavelengths are given by:

[TABLE]

where $\phi$ is now common across all wavelengths. In general, we refer to acquiring linear combinations of the array elements as coded aperture acquisition.

Yet another way of acquiring the samples in Figure 4 is to use coded frequency excitation where the excitation signal is a linear combination of various wavelengths. If the array elements are directly read out, then the measurements obtained can be interpreted as taking linear combinations of the samples along the radial lines in Figure 4.

A variant of aperture coding was considered in [9] where the authors propose using different subarrays at different times and then using interpolation techniques to acquire all the samples shown in Figure 4. Subsampling the array is equivalent to using binary codes on the aperture. However, their signal model is different from ours and hence the paper does not consider imaging with fewer than $M$ measurements.

We consider a particular signal model: we assume that the target to be imaged is range-limited. With this model, we study the number of array measurements required to image the scene when broadband excitation is used. We will show that coded apertures ((8) and (9)) can be used with broadband excitation to achieve highly efficient imaging of range-limited target scenes. Unlike standard methods that require $M$ direct read-outs or $M$ linear combinations of the array outputs, we will use far fewer “generic” linear combinations to image range-limited scenes. In other words, we show that when the target scene is range-limited, the coding matrices $\phi$ (or $\phi_{\lambda}$) can be highly underdetermined, without assuming any sparsity in the scene.

4 Signal model and degrees of freedom

In this section, we introduce our signal model and make some observations about the model. These observations, a signal acquisition scheme based on them, and theoretical guarantees on the proposed signal acquisition method are our main contributions.

4.1 Array measurements of range-limited scenes

Our signal model is that of a range-limited scene. A finite bandwidth and aperture allow us to observe only a part of the Fourier transform of the image. A finite range restricts the number of degrees of freedom of this observed region of the Fourier transform. We intend to take advantage of this to achieve a more efficient sampling of the bandlimited spectrum of the image and achieve hence faster imaging.

Briefly switching back to the setting of a continuous aperture, we can obtain an analog of (6) for range-limited targets. If $\mathcal{Q}$ represents the self-adjoint operator that truncates an image to a range limit $R$, the continuous aperture measurements can be modeled as

[TABLE]

The discrete array measurements can be modeled as

[TABLE]

We now demonstrate the effect of range-limitedness using a target scene with delta thickness present at a constant known range: where the scene has only one reflector per angle, with each reflector present at a constant known distance from the array. The underlying effect on the Fourier domain samples extends to scenes with a more general range limit.

Consider a scene with delta thickness at a constant range $R_{0}$ from the antenna array. The time domain outputs at different wavelengths can be modeled using (3) as:

[TABLE]

Considering just the amplitude as before, we have

[TABLE]

Now define

[TABLE]

The array outputs $y_{m}(\lambda)$ are then just ** samples of the same function** $g(\omega_{\tau})$ sampled uniformly in $[-D/\lambda,D/\lambda]$ (modulo known scaling factors). This is illustrated in Figure 5. In the 2D array case, the “slices” of the trapezoid corresponding to different excitation wavelengths sample a common function. As the range limit increases but remains finite, the functions sampled at different wavelengths start to differ, but still have limited degrees of freedom.

From the above discussion, it is clear that for scenes with delta thickness, collecting the full set of samples at the lowest excitation wavelength provides all available information. Collecting further samples at higher wavelengths offers no advantage. However, this redundancy can be used to collect fewer spatial samples but with broadband excitation. Similarly, collecting the full set of broadband measurements for scenes with higher but finite range limits results in a number of measurements greater than the number of degrees of freedom. It is thus natural to expect that the target scene can be reconstructed with a number of measurements $l\ll M$, owing to the limited number of degrees of freedom.

For computational purposes, we can discretize the target scene and the array imaging operator. For a scene with delta thickness at a depth $R_{0}$, let $x_{R_{0}}\in\mathbb{R}^{N}$ denote the target scene $x_{c}$ sampled uniformly with $N\geq M$. Similarly, the integral mapping the target scene to the array measurements at wavelength $\lambda$ can be discretized as a matrix $A_{R_{0},\lambda}\in\mathbb{C}^{M\times N}$ operating on $x_{R_{0}}$:

[TABLE]

With this notation in place, the array measurements can be expressed as

[TABLE]

Similarly, a more general scene with a range limit $R$ can be discretized as a vector of reflectivities as $x_{R}\in\mathbb{R}^{N\times D}$ where $N$ represents the number of discrete samples along the $\tau$ axis and $D$ the number of samples along the $r$ axis. Let the scene lie between the range limits $R_{\min}$ and $R_{\max}$. Define $d_{r}=\lfloor n/N+1\rfloor,\ n_{\tau}=\mod(n,N)$. Then the discretized array operator can be expressed as a matrix $A_{R}\in\mathbb{C}^{M\times ND}$ as

[TABLE]

and the array outputs can be expressed as

[TABLE]

Since our signal model considers only range-limited scenes, we drop the subscripts $R$ and $R_{0}$ for further discussion. We will denote antenna array measurements at excitation wavelength $\lambda$ as $y_{\lambda}\in\mathbb{C}^{M}$, the discretized target scene as vector $x_{0}\in\mathbb{R}^{ND}$ (or $\mathbb{R}^{N}$ for scenes with delta thickness) and the array imaging operator as $A_{\lambda}\in\mathbb{C}^{M\times ND}$ (or $\mathbb{C}^{M\times N}$ for scenes with delta thickness). $A_{\lambda}$ will incorporate knowledge of the target range profile, which is assumed to be known a priori. This helps us focus on the advantage of range-limitedness in itself. We later show how unknown range profiles can be handled algorithmically.

When the measurements at all the $K$ wavelengths are considered, we obtain the linear system

[TABLE]

The collection of array outputs at all the wavelengths lie in the column space of $A$. The effective dimension of this subspace determines the number of degrees of freedom in $y$. When the measurements are acquired with a coded aperture, we can model the measurements as

[TABLE]

Our main goal is to show that the target scene can be constructed using highly underdetermined matrices $\Phi$ with no loss in resolution compared to the reconstruction obtained using the full set of measurements $y$.

4.2 Degrees of freedom of range-limited scenes

In 4.1 we observed that bandlimited samples of range-limited target scenes have limited number of degrees of freedom. This is a consequence of simultaneous concentration of energy in both spatial and spectral domains. Signals with such a property can be well approximated by a number of basis functions that is proportional to the product of the area/volume of the spatial and spectral supports. This has been studied in a set of seminal papers by Slepian, Landau and Pollock [17, 18, 19, 29, 30, 20] and by Simons et.al., in [24, 25]. In essence, the range space of space-limiting, band-limiting operators such as that in (10) is approximately finite dimensional. This approximate dimension is the number of effective degrees of freedom of signals well concentrated in spatial and spectral domains. An efficient basis for the representation of such signals is the prolate basis [17]. Smaller the spatial and spectral supports, smaller is the number of degrees of freedom. See [21, 22, 23] for a quantitative non-asymptotic characterization of these properties in the discrete case.

Since the array imaging operator is an approximation of the continuous domain model in (10), its range space has a low dimensional structure. To illustrate this, we present the singular values of $A$ for target images with a finite range limit in Figure 6. The product of spatial and spectral supports in this case is approximately $2RD(\omega_{2}-\omega_{1})(\omega_{2}+\omega_{1})$. For these plots, an array with $213$ elements was used and samples were collected at $25$ wavelengths placed uniformly between $2$GHz and $4$GHz. The total number of samples collected is hence $5325$, but these samples lie in a subspace of dimension approximately only $640$ or $320$, for the two range limits considered. For the same antenna array, we show the singular values of the operators associated with target scenes that have delta thickness in Figure 7. As expected, such scenes lie subspaces of even smaller dimensions.

A special case of range-limitedness is when the target has delta thickness. Figure 7 shows the singular value decomposition of the array operator for two examples of such scenes: one where all the reflectors are at a constant known distance from the array center; the other where each reflector is at a known but different distance from the array center. We also show the Fourier transforms of the eigenvectors of $A^{T}A$ associated with the most and least significant eigenvalues in Figure 8.

The low dimensional structure in the array measurements forms the basis of our proposed imaging method. We show that this structure enables imaging with very few array measurements. In the following sections, we set up the reconstruction problem and provide theoretical guarantees for image reconstruction.

5 Coded aperture image reconstruction

5.1 Image reconstruction method

The limited number of degrees of freedom of broadband measurements leads us to ask the following question: what is the number of broadband measurements required to image range-limited target scenes? In order to answer that question, we first set up our reconstruction method. We do not assume any structure on the target scene such as sparsity, low total variation norm, or other structure apart from it being range-limited. Hence, with the measurement model used in (16), we solve the ordinary least squares problem

[TABLE]

Our standard of comparison is the ordinary least squares estimate obtained using the full data:

[TABLE]

The theoretical question we answer is: If each diagonal block $\phi$ in (16) is of size $l\times M$, what is the sample complexity $l$ that can achieve the same reconstruction result as (18) using (17)?

5.2 Sketched least squares systems

Problems of the type (17) are the subject of a recently popular area of research in the numerical linear algebra community, known as matrix sketching [5]. Sketching techniques use random projections of the rows (or columns) of a matrix $A$ to obtain approximations to the matrix. The new matrix obtained has much smaller dimensions and is used to solve the original linear algebra problem. Sketching can be applied to a multitude of problems such a linear regression, low rank matrix factorization, and subspace clustering among others.

If $y$ denotes measurements of an unknown signal $x$ observed through the linear operator $A$, then the least squares solution is given by (18). When $A$ is low rank ($\text{rank}(A)\ll\min(\tilde{M},N)$), $x$ can be solved for using fewer measurements. This forms the fundamental basis of the idea of sketching. Formally, if $S$ is a compressive matrix, sketching solves the following problem:

[TABLE]

where the subscript SLS stands for sketched least squares. The key point of course is that $SA$ is now a much smaller matrix and thus the solution can be computed quickly. This technique has generated a lot of interest in the context of large scale problems where the ambient dimensions are much bigger than the inherent dimensionality. In a sensing system, this translates to reconstructing a target signal from compressed measurements, as common in compressed sensing systems. However, we point out that unlike standard compressed sensing setup, we do not assume sparsity or compressibility in the signal $x$. Hence, we aim to solve (17) without imposing any structure on $x$.

Our proposed system can be modeled as a sketched linear system, but with a particularly structured sketching matrix:

[TABLE]

where $\Phi$ is a block diagonal matrix with $k$ repeated blocks

[TABLE]

We refer to such a matrix as a repeated block diagonal (RBD) matrix.

5.3 Linear algebraic interpretation of range-limitedness

The least squares program in (18) searches for an image that best explains the measurements obtained, and lies in the row space of $A$. The row space of $A$ is not only approximately low dimensional for range-limited images, but has a second tier of structure that makes aperture coding a very efficient way of obtaining fewer measurements. In particular, the relationship between the subspaces associated with the linear operators at various excitation wavelengths for range-limited images determines the number of aperture codes required.

Consider array imaging at a single wavelength $A_{\lambda_{i}}$. The least squares estimate in this case looks for a target in the row space of $A_{\lambda_{i}}$ that best explains the measurements. We denote this subspace as $\mathcal{S}_{i}$. Similarly, imaging using wideband excitation results in an estimate that lies in the union of subspaces $\mathcal{S}_{1},\mathcal{S}_{2},\cdots,\mathcal{S}_{k}$. The sample complexity of aperture coding depends highly on the relationship between these subspaces.

Let $A$ now denote a general $km\times n$ matrix of a finite rank $r$ such that $A=\begin{bmatrix}A_{1}^{T}&A_{2}^{T}&\cdots&A_{k}^{T}\end{bmatrix}^{T}$ where each $A_{i}$ is $m\times n$. Let the rank of $A_{i}$ be $r_{i}$. Without loss of generality, we can assume that $r_{i}\geq r_{j}$ for $i>j$, since the ordering of the row groups does not matter. We can then obtain the following factorization:

[TABLE]

where $C_{ij}\in\mathbb{R}^{m\times d_{j}}$, each $C_{ii}$ is full column rank when $d_{i}\neq 0$, and $V$ is an $r\times n$ orthonormal matrix. The factorization is such that the row space of $A_{1}$ is the span of the orthobasis $V_{1}$, the row space of $A_{2}$ is included in the span of $V_{1}$ and $V_{2}$. In general, $[V_{1}\ V_{2}\ \cdots V_{i}]$ includes an orthobasis for the row space of $A_{i}$. This factorization is equivalent to a block QR factorization of $A^{T}$ and can be obtained for any general matrix A. We are particularly interested in the $d_{i}$’s, as they capture the relationship between various subspaces. The diagonal blocks $C_{ii}$’s represent the energy in the subspace orthogonal to the union of the row spaces of the previous blocks $A_{1},\cdots,A_{i-1}$. Hence, smaller values of $d_{i}$ indicate that the subspaces have a high degree of overlap.

A special case of high overlap is when the row spaces have a nested structure: $\operatorname{row}(A_{i-1})\subset\operatorname{row}(A_{i})\ \forall i=2,\cdots,K$. In this case, the off-diagonal blocks $C_{ij},i<j$ and the orthogonal blocks $V_{1},\cdots,V_{i-1}$ capture a significant part of $\operatorname{row}(A_{i})$. For low rank systems, this naturally leads to smaller values of $d_{i}$. In contrast, when the row spaces are all almost orthogonal, the $d_{i}$’s are all large and the off-diagonal blocks $C_{ij},i\neq j\approx 0$. We will later show that the broadband array imaging operator has the nested subspace structure for certain range-limited scenes.

Let us now relate the above factorization to the context of imaging. To begin, suppose that we use wavelengths up to $\lambda_{i}$. Then, $d_{i+1}$ represents the rank of the update required to incorporate information from a new, lower wavelength $\lambda_{i+1}$. It represents the innovation added by the measurements at the new wavelength.

As the range limit decreases, the overlap between the subspaces $\mathcal{S}_{i}$ increases, increasing the redundancy across wavelengths. This leads to the possibility of higher subsampling rates in the physical array domain. In the limiting case of only one reflector per angle, the subspaces have a nested structure. As we will observe in Section 6, this plays a crucial role in determining the number of aperture codes required for successful imaging.

Our theoretical results formally state the effect of the relationship between the subspaces $S_{i}$ on the number of aperture codes required for imaging. Theorem 6.1 provides a non-trivial estimate of the number of measurements needed for a given excitation bandwidth. Theorem 6.4 provides conditions under which a given set of excitation wavelengths $\lambda_{1},\cdots,\lambda_{K}$ allow the number of spatial measurements to be reduced by a factor of $k$. Using these conditions and a given bandwidth of excitation, one can choose the set of excitation wavelengths and very few coded measurements to achieve imaging with no loss in resolution.

6 Signal recovery from aperture coded measurements

In the previous section, we set up the aperture coding problem as a sketched least squares problem that has a particular structure dictated by the physical problem of array imaging. In this section, we derive mathematical guarantees for such sketched systems and provide estimates of the required sample complexity.

We start by reviewing the conditions that any general sketching operator has to satisfy in order for the solution to the sketched least squares problem $x_{\text{SLS}}$ to be close to the solution of the original ordinary least squares solution $x_{\text{LS}}$. In the noiseless case

[TABLE]

where $A=U\Sigma V^{*}$ is the SVD of the linear operator and $A^{\dagger}$ denotes the pseudoinverse of $A$. (23) shows that the least squares solution is a projection of the true solution onto the row space of $A$. Hence, any sketching operator $\Phi$ should preserve the row space of $A$. It has been well established in literature that a number of random projections of the rows of $A$ greater than or equal to its rank capture the row space in case of exactly low rank matrices [27]. When $A$ of size $m\times n$ and $\text{rank}(A)=r\ll\min(m,n)$, if $\Phi$ is a $l\times m$ dense standard normal random matrix with $l\geq r$, then $\operatorname{row}(A)\subset\operatorname{row}(\Phi A)$. Since by construction we also have $\operatorname{row}(\Phi A)\subset\operatorname{row}(A)$, we have

[TABLE]

Let

[TABLE]

The least squares estimate from the sketched measurements $\phi y$ is same as that from the full observation $y$, since

[TABLE]

This idea forms the basis of using sketched measurements to solve a least squares problem. For any sketching matrix $\phi$, a necessary and sufficient condition in the noiseless case is $\|(I-P_{(\phi A)^{*}})A^{*}\|=0$.

Our goal is to replicate the result that $l\geq r$ is sufficient, but for an RBD matrix. For such a matrix (shown in (21)), the equivalent result would be that a total number of measurements $Kl\geq r$ suffice to capture the row space of the matrix $A$. However, due to the highly structured nature of a block diagonal matrix, such a result does not hold uniformly for all matrices $A$. We analyze the conditions on $A$ under which such a result holds and show that array imaging matrices do obey these conditions, thus allowing for spatial subsampling.

6.1 Least squares with a block diagonal sketching matrix: General case

RBD matrices obtain localized random projections: they take linear combinations of only a subset of the rows. In this section, we provide guarantees on the the error $\|(I-P_{(\Phi A)^{*}})A^{*}\|$ when $\Phi$ is an RBD matrix. The focus will be on the sample complexity $l$ required to drive this error to [math] with high probability.

It is immediately clear how to achieve this when we take $\displaystyle l\geq\max_{i}\text{rank}(A_{i})$. Let $A=[A_{1}^{T}\ A_{2}^{T}\cdots\ A_{K}^{T}]^{T}$ where each $A_{i}$ is of size $m\times n$ and has rank $r_{i}$. Let $Y=\Phi A$. Hence

[TABLE]

Since

[TABLE]

and $\|(I-P_{Y_{i}^{*}})A_{i}^{*}\|=0$ for $l\geq r_{i}$, we obbtain

[TABLE]

for $l\geq\max_{i}r_{i}$.

Compared to using a dense random matrix, this can be worse by a factor of $K$. Intuitively, this straightforward application of results from [27] leads to capturing of the subspace spanned by each group of rows $A_{i}$ individually, without considering the overlap between the subspaces. This is addressed in our first analytical result (Theorem 6.1), which provides a simple but non-trivial estimate of the number of random projections required. We then improve this result in Theorems 6.3 and 6.4.

Theorem 6.1**.**

*For a given matrix $A$ of size $Km\times n$, let the $d_{i}$ be defined as in (22). Let $\Phi$ be a block diagonal matrix with repeated diagonal block $\phi$ of size $l\times m$ and whose entries are chosen i.i.d. from the standard normal distribution. Let $Y=\Phi A$. Define $d_{0}=\max_{i}d_{i}$. For $l\geq d_{0}$, $\|(I-P_{Y^{*}})A^{*}\|=0$ with probability $1$. *

Proof 6.2**.**

We have

[TABLE]

*Then, if each diagonal block is full rank, the matrix $\Phi C$ and hence $\Phi CV^{T}$ is full rank, since $V$ is just an orthonormal matrix. Consider each diagonal block of $\Phi C$, $\phi C_{ii}$. Since $l\geq d_{0}$, $\operatorname{rank}(\phi C_{ii})=d_{i}$ with probability $1$. Since the rank of any block triangular matrix is at least the sum of ranks of the diagonal blocks, $\operatorname{rank}(\Phi C)=\operatorname{rank}(\Phi CV^{T})=\sum_{i}d_{i}=r$. Since $\operatorname{row}(\Phi CV^{T})\subset\operatorname{row}(CV^{T})$ and $\operatorname{rank}(\Phi CV^{T})=\operatorname{rank}(CV^{T})$, $\operatorname{row}(\Phi CV^{T})=\operatorname{row}(CV^{T})$. Hence, $P_{Y^{*}}=P_{A^{*}}$ and the conclusion follows. *

To understand this theorem, we can first consider the intuition behind random projections: owing to the randomness of the projections, a linearly independent set of vectors in the row/column space of the original matrix is obtained and hence the subspace is captured. With an RBD sketching matrix, one can only obtain random projections within the subspace of each group of rows. However, if the subspaces spanned by these row groups overlap, we may obtain a linearly independent set of vectors that capture the row space of the whole matrix by using a few random projections of each row group. The factorization in (22) captures this dependence between the subspaces spanned by the row groups.

The bound in Theorem 6.1 can be tight or loose depending on the candidate matrix $A$. The bound is tight when the row spaces span non-overlapping subspaces, since in this case, random projections within a row group do not provide any information about the other row groups. When the subspaces do overlap the bound can be significantly improved, even though it is already better than the trivial bound $l\geq\max_{i}r_{i}$. For example, consider a scenario where the subspaces spanned by the $A_{i}$’s are nested in the subspaces spanned by the $A_{j}$’s for all $i<j$. In such a case, by obtaining random projections within the span of $A_{i}$, we are also already obtaining random projections in the subspace spanned by $A_{j}$ for all $i<j$. This allows for $l$ to be much smaller, even if some of the $d_{i}$ are large, contrary to what Theorem 6.1 predicts. As we note in the next subsection, this is precisely the case for the linear operator associated with range-limited images.

Although our result concerns exactly low-rank matrices, it can be generalized to matrices that are approximately low rank using perturbation theory for projection matrices [31]. In general, for matrices with numerical ranks much smaller than the ambient dimensions, the subspaces spanned by the $A_{i}$’s will overlap. This overlap reduces the number of random vectors needed per block.

Figure 9 shows the error $\|(I-P_{Y^{*}})A^{*}\|$ and the spectra of the matrices $A$ and $\Phi A$ for a randomly generated test matrix of size $2000\times 1000$, with $\max_{i}(r_{i})\approx 300$ and $d_{0}=110$. We observe empirically that the row space was captured with $l\approx 40$.

6.2 Least squares with a block diagonal sketching matrix: Overlapping subspaces

Consider a 1D array and a 2D image. For an image with delta thickness with reflectors at a known fixed depth $R_{0}$, the rows of $A_{\lambda_{i}}$ are unit-length sinusoids in the frequency range $[-D/\lambda_{i},\ D/\lambda_{i}]$, with a modulation term $e^{-j2\pi R_{0}/\lambda_{i}}$. The row space of $A_{i}$ is therefore well approximated by the first $2D/\lambda_{i}$ discrete prolate spheroidal sequences (DPSS) [18]. As the wavelength progresses from the highest $\lambda_{\max}$ to the lowest $\lambda_{\min}$, the subspace spanned by the rows of $A_{\lambda_{i}}$ is nested in the subspace spanned by $A_{\lambda_{j}}$ for all $j>i$. This is shown in Figure 10 where the full operator with all wavelengths has a row space of the same dimension as the operator at only the highest frequency.

Although Theorem 6.1 is appealing because of its simplicity, it can lead to suboptimal estimates of $l$. For example, if $\lambda_{\min}=7.5$cm and $\lambda_{\max}=15$cm, one could only hope for a reduction in the number of measurements by a factor of $2$. But as we have already seen in Figure 1, we can easily achieve better subsampling.

The following two theorems provide much stronger results: they directly address the question of when $l\geq r/K$ random projections suffice to achieve $\|(I-P_{Y^{T}})A^{T}\|=0$. In the context of array imaging, they address the question of when using $K$ different excitation wavelengths can offer the luxury of imaging with roughly only $M/K$ measurements. Their proofs are deferred to the appendix.

Theorem 6.3**.**

*Let A=\bigl{[}\begin{smallmatrix}A_{1}\\ A_{2}\end{smallmatrix}\bigr{]}, with $\operatorname{row}(A_{1})\subseteq\operatorname{row}(A_{2})$. Assume $A_{2}$ is full row rank and define $U=A_{1}A_{2}^{T}$. If the entries of $\phi$ are drawn from a continuous distribution, with probability 1, M=\bigl{[}\begin{smallmatrix}\phi A_{1}\\ \phi A_{2}\end{smallmatrix}\bigr{]} is full row rank for any $l\leq n/2$ if and only if no real eigenvalue of $U$ has an algebraic multiplicity greater than $n/2$. Consequently, for $l\geq n/2$, $\|(I-P_{M^{T}})A^{T}=0\|$ with probability 1. *

Theorem 6.3 provides the necessary and sufficient conditions on an ensemble of two matrices with a nested subspace structure under which only $r/2$ projections of each block are sufficient to capture the row space. The condition prohibits the existence of an invariant subspace of dimension greater the $n/2$ that is common to both $A_{1}$ and $A_{2}$. Otherwise, $\phi A_{1}$ and $\phi A_{2}$ will both have a component along this subspace for some $l<n/2$, resulting in a loss of linear independence.

The eigenvalue distribution for the array imaging matrices is shown in Figure LABEL:sub@fig:eigU when $\lambda_{\min}=7.5$cm and $\lambda_{\max}=10$cm. It is clear that the matrices meet the required condition and hence only $M/2$ measurements suffice, in contrast to what is predicted by Theorem 6.1, which would be at least $3M/4$.

In our next result, we extend the result to the case with more than two blocks and provide a sufficient condition on the ensemble of $K$ matrices $\{A_{i}\}$ under which only $l\geq r/K$ random projections per block can capture the full row space. To do this, we define the following matrices:

[TABLE]

where $V_{\mathcal{S}}$ is any orthonormal matrix of size $l\times m$.

Theorem 6.4**.**

*Given an ensemble of $K$ matrices $\{A_{i}\}$ for $i=1,\cdots,K$, each of size $m\times n$ and a matrix $\phi\in\mathbb{R}^{l\times m}$ with entries drawn from the standard normal distribution, $M$ as defined above is full row rank if there exists an orthonormal basis $V_{\mathcal{S}}\in\mathbb{R}^{l\times m}$ such that the $kl\times n$ size matrix $\widehat{M}$ has full row rank. Consequently, for $l=r/K$, $\|(I-P_{M^{T}})A^{T}\|=0$ with probability 1. *

Intuitively, Theorem 6.4 requires that there is at least one subspace of dimension $l=r/K$, which when projected onto the matrices $A_{i}$ results in a set of $K$ linearly independent subspaces. In Figure LABEL:sub@fig:hist, we show the histogram of the smallest singular values of the matrix $\widehat{M}$ for the array imaging operator with $K=8$ excitation wavelengths placed uniformly between $\lambda_{\min}=7.5$cm and $\lambda_{\max}=15$cm, for $1000$ realizations of randomly chosen orthonormal basis $V_{\mathcal{S}}$. In this case, the number of array elements was $M=213$ and the scene considered had delta thickness. Hence, approximately only $30$ spatial measurements suffice in imaging any such scene. As the range extent of the images increases, the nested structure in the row spaces ceases to exist. However, the subspaces are still have a high degree of overlap and a number of aperture codes much smaller than the number of conventional beams used suffice for imaging.

The conditions stated in Theorems 6.3 and 6.4 ensure that the lack of diversity among the diagonal blocks of the RBD matrix is compensated for by the data matrix $A$ itself. As we show in further sections, the array imaging operator satisfies these conditions, thus lending itself to spatial subsampling.

7 Extensions to aperture coded imaging

The proposed aperture coded imaging system requires obtaining linear combinations of the array outputs at different excitation wavelengths. The analysis presented in the previous section assumed that the weights for these linear combinations are drawn from a continuous distribution such as the standard normal distribution and also assumed noise-free conditions. In this section, we wish to point out an acquisition method that follows a similar spirit but can be practically more appealing and we also model noisy imaging settings. Although we do not provide analysis for these cases, our simulations, shown in the next section, demonstrate that the fundamental idea of trading spatial measurements with excitation bandwidth for range limited images can still offer significant gains.

7.1 Imaging with subsampled array

The observation that the row spaces of the operator at different wavelengths overlap to a high degree for range-limited images suggests that we can also directly subsample the array and obtain direct measurments without any loss in resolution. Subsampled arrays can be thought of as aperture coded arrays with binary codes. Again drawing a parallel to the literature in randomized numerical linear algebra, subsampling the array is similar to sampling a few rows and columns of a matrix to provide an approximation. (See [32] and the references therein.) A standard approach in the linear algebraic community is to use randomized subsampling. A similar scheme can be applied to antenna arrays: sample only a fraction of randomly chosen array elements with broadband excitation. This results in a set of array measurements at different wavelengths that share the same subsampling pattern.

There could also be more principled subsampling approaches for the specific application of broadband array imaging. Considering Figure 5, one might expect that sampling a few terminal elements and some elements at the center may be sufficient to obtain enough samples to reconstruct the common function. Providing theoretical justification for how many array elements are required could be an interesting direction for future research. In Figure 12, we show the binary mask resulting from such an approach described for a 2D array. In the case of a scene with delta thickness present at a known constant range, measurements at each wavelength only samples the innovation that the wavelength adds compared to other larger wavelengths.

7.2 Imaging in a noisy scenario

The theorems from previous section address the case of noiseless imaging. However, we show in our experiments that aperture coding is also robust to noise. Stable least squares reconstruction is commonly achieved using Tikhonov regularization. Noisy measurements are modeled as

[TABLE]

where $e$ represents the noise vector. $x$ can be estimated by solving

[TABLE]

where $\delta$ is a regularization parameter. Robust reconstruction can be obtained using the right choice of the regularization parameter. For aperture coding, we assume the following model for noisy measurements, just as in traditional beamforming methods:

[TABLE]

Again, we estimate $x$ by solving

[TABLE]

The choice of $\delta$ and $\hat{\delta}$ affect the reconstruction performance. While increasing the regularization parameter reduces the effect of noise on the reconstruction, it also leads to an increase in the signal reconstruction performance. The choice of the parameter hence leads to a trade-off. In our simulations, we show that for a given choice of $\delta$, we can choose $\hat{\delta}$ to match the signal reconstruction performance and the noise performance of traditional imaging.

We use the following definitions to quantify the performance of aperture coding in a noisy scenario:

[TABLE]

[TABLE]

where $(S)^{\dagger}_{\operatorname{Tikh}}=(S^{T}S+\delta I)^{-1}S^{T}$ denotes the Tikhonov regularized pseudo-inverse of the matrix $S$.

8 Experiments

We now provide experimental results to show the effectiveness of aperture coding. Various experiments were conducted: aerture coding simulations were conducted with images at a constant depth from the antenna array, for flat images parallel to the 2D array and delta thickness images that are multi-depth with known and unknown range profiles. The next set of experiments deal with subsampling the array for the same class of images. These are followed by simulations in the presence of noise. We then provide simulation results for objects whose range limits are not thin and span a range of $15$cm. The following array parameters were used in all the simulations: an array of $40\times 40$ elements were used with 15 excitation wavelengths placed regularly in the bandwidth of $7.5$cm and $15$cm. The elements were placed at half the smallest wavelength. The scene was assumed to be within the angular span of $[-\pi/4,\ \pi/4]$ in both directions. For this configuration the number of beams/ measurements needed at each frequency for conventional imaging is about 1100. Quantitative error values for all the experiments are given in Table 8.

A. Constant range and multi-depth images

Figure 1 shows reconstruction results for images at a constant, known range. Conventional beamforming requires about 1100 beams in this case to scan over the entire image. It is clear that similar reconstruction performance can be obtained using as few as 80 beams with wideband excitation. The relative reconstruction error is almost negligible in each case (shown in Table 8).

Aperture coding is also effective for multi-depth images. A scene with three segments, each at a different depth was considered. The depth map of the scene used in the simulations is shown in Figure 13. In Figure 14, we present simulation results when the depth profile was assumed to be known a priori. Again only about 80 beams are sufficient to get good quality reconstructions.

We further explore the performance of coded aperture imaging when the depth profile is unknown. Boufounos in [33] describes a method to infer the depth profile of an image using CoSaMP algorithm. The algorithm uses full array measurements in the least squares step. We replace this step with sketched least squares and coded measurements. Simulation results are shown in Figure 15. We note that aperture coding still performs well, with the relative error being negligible even for 160 generic beams.