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 . Nevertheless, applications such as autonomous vehicles, depth sensing, gesture recognition  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.
We use the standard linear model for array imaging in a setting with a single transmitter and multiple receiver elements. We use to denote the number of antenna array elements and 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 , , can be modeled as
where denotes the sampled target scene and is the linear operator ( matrix) mapping the scene to the measurements. Collecting broadband measurements is equivalent to collecting for . We note that all can be obtained simultaneously (using only
“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:
where is typically an invertible matrix. In traditional beamforming the weights are chosen to induce spatial selectivity where each of the 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 ofcompressed 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 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  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
where is a highly structured compresssive matrix. We show that it provides a solution equivalent to that of solving the higher dimensional problem
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 sensors with sensors placed cm apart and a target scene having a span of 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 beams at wavelength cm. By introducing bandwidth in the excitation signal, we show that the scene can be imaged with as few as spatial linear combinations. This is illustrated in Figure 1. For target scenes with higher range limits, standard imaging methods also need bandwidth for imaging . 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 . 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,  proposes to increase the number of resolvable directions to by using an array with 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,  imposes sparsity based regularization to solve the ill-posed radio-interferometric imaging problem. A similar theme is also followed in  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 , 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  where the Fourier basis is used as a sketching matrix for dimensionality reduction in interferometry. But the sketching matrix in  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 placed on the -axis, spanning . 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 . The system consists of only one transmitting element which is co-located with the receiver at the array center. Let be the continuous time index. For an excitation signal , the signal received at the array location is for a unit-strength reflector at . For a general reflectivity map, the narrowband response at this location for the excitation signal is
By making the substitution , the complex amplitude of the signal received at location for excitation wavelength can be written as
where , and denotes the Fourier transform of . This shows that the antenna aperture measures the Fourier transform of (a warped version of) the target scene. can be obtained by computing the Fourier transform of the temporal signal 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 , the accessible interval in the Fourier domain for is limited to .
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 , where is the roundtrip distance to the array center, is the angle with respect to the Y-Z plane, and is the angle with respect to the X-Z plane, as shown in Figure (b)b. The scene reflectivity is denoted as where and . At excitation wavelength , the 2D array outputs are samples of the Fourier transform of the scene sampled in the region bounded by , where and 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 . If the broadband signal used is , then the received signal at location is
The complex amplitudes at different wavelengths, can then be obtained as the Fourier coefficients of 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 represents the continuous domain Fourier transform and represents the bandlimiting operator supported only in this accessible window, the continuous aperture measurements can be modeled as
where 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 discrete antenna elements placed uniformly at coordinates . For ease of explanation, we consider a discrete set of excitation wavelengths . 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 where ,
can be denoted by a vector. The set of measurements at each wavelength is denoted by .
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 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 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 . Instead of acquiring direct read-outs, 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 , the measurements made in time domain are
for . The vector of complex amplitudes is then given by
where is in general an well-conditioned matrix whose element is . 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:
where 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 
where the authors propose using different subarrays at different times and then using interpolation techniques to acquire all the samples shown in Figure4. 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 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 direct read-outs or 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 (or ) 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 represents the self-adjoint operator that truncates an image to a range limit , the continuous aperture measurements can be modeled as
The discrete array measurements can be modeled as
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 from the antenna array. The time domain outputs at different wavelengths can be modeled using (3) as:
Considering just the amplitude as before, we have
The array outputs are then just samples of the same function sampled uniformly in (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 , 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 , let denote the target scene sampled uniformly with . Similarly, the integral mapping the target scene to the array measurements at wavelength can be discretized as a matrix operating on :
With this notation in place, the array measurements can be expressed as
Similarly, a more general scene with a range limit can be discretized as a vector of reflectivities as where represents the number of discrete samples along the axis and the number of samples along the axis. Let the scene lie between the range limits and . Define . Then the discretized array operator can be expressed as a matrix as
and the array outputs can be expressed as
Since our signal model considers only range-limited scenes, we drop the subscripts and for further discussion. We will denote antenna array measurements at excitation wavelength as , the discretized target scene as vector (or for scenes with delta thickness) and the array imaging operator as (or for scenes with delta thickness). 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 wavelengths are considered, we obtain the linear system
The collection of array outputs at all the wavelengths lie in the column space of . The effective dimension of this subspace determines the number of degrees of freedom in . When the measurements are acquired with a coded aperture, we can model the measurements as
Our main goal is to show that the target scene can be constructed using highly underdetermined matrices with no loss in resolution compared to the reconstruction obtained using the full set of measurements .
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 . 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 offor target images with a finite range limit in Figure 6. The product of spatial and spectral supports in this case is approximately . For these plots, an array with elements was used and samples were collected at wavelengths placed uniformly between GHz and GHz. The total number of samples collected is hence , but these samples lie in a subspace of dimension approximately only or , 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 ofassociated 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
Our standard of comparison is the ordinary least squares estimate obtained using the full data:
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 . Sketching techniques use random projections of the rows (or columns) of a matrix
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 denotes measurements of an unknown signal observed through the linear operator , then the least squares solution is given by (18). When is low rank (), can be solved for using fewer measurements. This forms the fundamental basis of the idea of sketching. Formally, if is a compressive matrix, sketching solves the following problem:
where the subscript SLS stands for sketched least squares. The key point of course is that 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 . Hence, we aim to solve (17) without imposing any structure on .
Our proposed system can be modeled as a sketched linear system, but with a particularly structured sketching matrix:
where is a block diagonal matrix with repeated blocks
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 . The row space of 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 . The least squares estimate in this case looks for a target in the row space of that best explains the measurements. We denote this subspace as . Similarly, imaging using wideband excitation results in an estimate that lies in the union of subspaces . The sample complexity of aperture coding depends highly on the relationship between these subspaces.
Let now denote a general matrix of a finite rank such that where each is . Let the rank of be . Without loss of generality, we can assume that for , since the ordering of the row groups does not matter. We can then obtain the following factorization:
where , each is full column rank when , and is an orthonormal matrix. The factorization is such that the row space of is the span of the orthobasis , the row space of is included in the span of and . In general, includes an orthobasis for the row space of . This factorization is equivalent to a block QR factorization of and can be obtained for any general matrix A. We are particularly interested in the ’s, as they capture the relationship between various subspaces. The diagonal blocks ’s represent the energy in the subspace orthogonal to the union of the row spaces of the previous blocks . Hence, smaller values of 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: . In this case, the off-diagonal blocks and the orthogonal blocks capture a significant part of . For low rank systems, this naturally leads to smaller values of . In contrast, when the row spaces are all almost orthogonal, the ’s are all large and the off-diagonal blocks . 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 . Then, represents the rank of the update required to incorporate information from a new, lower wavelength . It represents the innovation added by the measurements at the new wavelength.
As the range limit decreases, the overlap between the subspaces 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 on the number of aperture codes required for imaging. Theorem 1 provides a non-trivial estimate of the number of measurements needed for a given excitation bandwidth. Theorem 3 provides conditions under which a given set of excitation wavelengths allow the number of spatial measurements to be reduced by a factor of . 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 to be close to the solution of the original ordinary least squares solution . In the noiseless case
where is the SVD of the linear operator and denotes the pseudoinverse of . (23) shows that the least squares solution is a projection of the true solution onto the row space of . Hence, any sketching operator should preserve the row space of . It has been well established in literature that a number of random projections of the rows of greater than or equal to its rank capture the row space in case of exactly low rank matrices . When of size and , if is a dense standard normal random matrix with , then . Since by construction we also have , we have
The least squares estimate from the sketched measurements is same as that from the full observation , since
This idea forms the basis of using sketched measurements to solve a least squares problem. For any sketching matrix , a necessary and sufficient condition in the noiseless case is .
Our goal is to replicate the result that 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 suffice to capture the row space of the matrix . However, due to the highly structured nature of a block diagonal matrix, such a result does not hold uniformly for all matrices . We analyze the conditions on 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 when is an RBD matrix. The focus will be on the sample complexity required to drive this error to
with high probability.
It is immediately clear how to achieve this when we take . Let where each is of size and has rank . Let . Hence
and for , we obbtain
Compared to using a dense random matrix, this can be worse by a factor of . Intuitively, this straightforward application of results from  leads to capturing of the subspace spanned by each group of rows individually, without considering the overlap between the subspaces. This is addressed in our first analytical result (Theorem 1), which provides a simple but non-trivial estimate of the number of random projections required. We then improve this result in Theorems 2 and 3.
Then, if each diagonal block is full rank, the matrix and hence is full rank, since is just an orthonormal matrix. Consider each diagonal block of , . Since , with probability . Since the rank of any block triangular matrix is at least the sum of ranks of the diagonal blocks, . Since and , . Hence, 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 1 can be tight or loose depending on the candidate matrix . 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 . For example, consider a scenario where the subspaces spanned by the ’s are nested in the subspaces spanned by the ’s for all . In such a case, by obtaining random projections within the span of , we are also already obtaining random projections in the subspace spanned by for all . This allows for to be much smaller, even if some of the are large, contrary to what Theorem 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 . In general, for matrices with numerical ranks much smaller than the ambient dimensions, the subspaces spanned by the ’s will overlap. This overlap reduces the number of random vectors needed per block.
Figure 9 shows the error and the spectra of the matrices and for a randomly generated test matrix of size , with and . We observe empirically that the row space was captured with .
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 , the rows of are unit-length sinusoids in the frequency range , with a modulation term . The row space of is therefore well approximated by the first discrete prolate spheroidal sequences (DPSS) . As the wavelength progresses from the highest to the lowest , the subspace spanned by the rows of is nested in the subspace spanned by for all . 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 1 is appealing because of its simplicity, it can lead to suboptimal estimates of . For example, if cm and cm, one could only hope for a reduction in the number of measurements by a factor of . 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 random projections suffice to achieve . In the context of array imaging, they address the question of when using different excitation wavelengths can offer the luxury of imaging with roughly only measurements. Their proofs are deferred to the appendix.
Let , with . Assume is full row rank and define . If the entries of are drawn from a continuous distribution, with probability 1, is full row rank for any if and only if no real eigenvalue of has an algebraic multiplicity greater than . Consequently, for , with probability 1.
Theorem 2 provides the necessary and sufficient conditions on an ensemble of two matrices with a nested subspace structure under which only projections of each block are sufficient to capture the row space. The condition prohibits the existence of an invariant subspace of dimension greater the that is common to both and . Otherwise, and will both have a component along this subspace for some , resulting in a loss of linear independence.
The eigenvalue distribution for the array imaging matrices is shown in Figure LABEL:sub@fig:eigU when cm and cm. It is clear that the matrices meet the required condition and hence only measurements suffice, in contrast to what is predicted by Theorem 1, which would be at least .
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 matrices under which only random projections per block can capture the full row space. To do this, we define the following matrices:
where is any orthonormal matrix of size .
Given an ensemble of matrices for , each of size and a matrix with entries drawn from the standard normal distribution, as defined above is full row rank if there exists an orthonormal basis such that the size matrix has full row rank. Consequently, for , with probability 1.
Intuitively, Theorem 3 requires that there is at least one subspace of dimension , which when projected onto the matrices results in a set of linearly independent subspaces. In Figure LABEL:sub@fig:hist, we show the histogram of the smallest singular values of the matrix for the array imaging operator with excitation wavelengths placed uniformly between cm and cm, for realizations of randomly chosen orthonormal basis . In this case, the number of array elements was and the scene considered had delta thickness. Hence, approximately only 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.
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  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
where represents the noise vector. can be estimated by solving
where 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:
Again, we estimate by solving
The choice of and 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 , we can choose 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:
where denotes the Tikhonov regularized pseudo-inverse of the matrix .
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 cm. The following array parameters were used in all the simulations: an array of elements were used with 15 excitation wavelengths placed regularly in the bandwidth of cm and cm. The elements were placed at half the smallest wavelength. The scene was assumed to be within the angular span of 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 1.
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 1).
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  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.
|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|
C. Imaging with noise
As claimed before, aperture coded imaging can still perform well in the presence of noise. We simulated an imaging scenario at dB input SNR level. In Figure 16, we show reconstruction results of these simulations. The corresponding signal reconstruction error values and output SNR are given in Table 1. It is clear that aperture coded measurements do not result in any degradation in the presence of noise. This shows that the sketched matrix is not only full rank, but also has stable eigenvalues.
D. Imaging scenes with higher range limits
Aperture coding can be highly effective even when imaging scenes with higher range limits. As predicted by our theorems, when the range limit increases, more measurements are required, but it can still less be than that used in conventional imaging. We demonstrate this in our next set of simulations. We consider an object made up of five layers as shown in Figure (a)a. These layers lie in a region of width cm, in the far-field of a 2D array at a depth of m. We again use 15 excitation wavelengths between cm and cm. In Figure 17, we present reconstruction results by using 640, 480, and 320 aperture codes in place of 1100 beams used in conventional imaging.
In this paper, we identify that broadband array imaging of range-limited target scenes is a particular case of spatio-spectral concentration, which results in the reconstructed image having a limited number of degrees of freedom. However, any image acquisition method that takes advantage of this low dimensional structure has to also be practically feasible. We show both theoretically and using simulations that existing array architectures such as aperture codes and subsampled arrays offer a way to make use of the low dimensionality to perform array imaging with far fewer spatial measurements than conventional methods when broadband excitation is used.
To establish the theoretical guarantees, we modeled the array imaging problem as a sketched least squares problem with a particular sketching matrix that has a block diagonal structure. We show that even with such structured randomness, the sketching matrix can be as small/compressive as a more generic dense sketching matrix. In terms of array imaging, this implies that existing array architectures can be used along with broadband excitation to image range-limited target scenes with very few spatial measurements.
Appendix A Proof of Theorem 2
The main tool we use to prove Theorems 2 and 3 is polynomial identity testing: if a polynomial of a finite total degree is said to be identically zero if the coefficient of every monomial term is zero. Any polynomial which is not identically zero, when evaluated at a random point drawn from a continuous distribution, is non-zero with probability 1, since the set of roots of the polynomial has measure zero with respect to the field of real numbers. In conclusion, a multivariate polynomial evaluated at a multi-dimensional point drawn from a continuous distribution is non-zero with probability 1, if any of the polynomial coefficients are non-zero.
To show that a given matrix is of rank , we choose a suitable submatrix of size and show that it has a determinant with at least one non-zero coefficient and therefore not identically zero. For random projections whose entries are drawn from the standard normal distribution, we hence need to show that the sketched matrices have submatrices of desirable sizes that have a non-zero determinant with probability 1.
For ease of notation, we prove the result on the transposed matrices: , , , and is the sketching matrix with entries drawn from the standard normal distribution. For ease of explanation, we assume that and are orthonormal matrices: and . We later explain how our proof holds for any two general matrices.
Consider the matrix where
is an orthogonal matrix. Ifis full column rank, then is full column rank. Any orthogonal matrix can be decomposed as
where is an orthobasis and is a block diagonal matrix with or blocks. The blocks are of the form
with being a pair of complex conjugate eigenvalues of , . The diagonal blocks are equal to and are also eigenvalues of . (32) is referred to as the canonical decomposition of .
Notice that . Since is orthogonal, is also a matrix with i.i.d. standard normal variables due to the rotational invariance of the standard normal distribution. Thus we can directly work with the matrix . The following lemma provides the necessary condition on the multiplicity of any real eigenvalue for to be full column rank.
Let there be a real eigenvalue of that has an algebraic multiplicity . Then for some , is not full column rank.
Let be a real eigenvalue with algebraic multiplicity and let . Then
where is the submatrix formed by the first rows of and is the submatrix formed by the last rows and columns of . Let , . Then, the submatrix has a null space of dimension and the submatrix has an dimensional null space given by the range of . Since , . The result follows.
We now prove some lemmas that provide sufficient conditions on for which is full column rank. We state the lemmas along with their proofs and then use them to prove Theorem 2.
Let there be pairs of complex conjugate eigenvalues of with non-zero imaginary part. Then for any , has full column rank.
Let be expressed as
where . Let denote the th element of . Rearranging the columns of , we obtain
Expanding the determinant of , the coefficient of the term is . For any , has a term of the form with coefficient