DOA Estimation in Antenna Arrays: Super-Resolution Methods in a Common Space

1. Introduction

Direction of Arrival (DOA) estimation is a foundational technique in array signal processing, with critical applications in radar, sonar, and modern wireless communications [1]. To achieve higher angular resolution beyond traditional limits, researchers have introduced various super-resolution (sub-Rayleigh) algorithms. Prominent examples include Maximum Entropy (ME) [2], Capon [3], Pisarenko [4], MUSIC [5], ESPRIT [6], and VESPA [6]. These methods usually fall into two views: linear prediction or signal-noise subspaces. However, strong links exist within and between these groups. This paper presents a common view of angular super-resolution methods as linear prediction of signals beyond an antenna aperture. This framework shows that shifting between methods simply involves changing one or more of these three parameters. It also allows us to derive new estimates that lie anywhere between existing methods.

This approach enables universal algorithms to produce spatial spectrum estimates across multiple methods. We demonstrate this by using the J orthogonalization algorithm to obtain estimates of the Pisarenko, maximum entropy, Capon and MUSIC spectra, as well as any intermediate estimates positioned between these methods.

The paper is organized as follows. Section 2 examines angular super-resolution methods using linear signal prediction beyond an antenna aperture, ranging from maximum entropy (ME) to minimum entropy (Pisarenko).

Section 3 examines multichannel linear prediction using subarrays with varying entropy; within this framework, the minimum and maximum entropy approaches correspond to the original and approximate ESPRIT methods, respectively.

In Section 4 we demonstrate that VESPA-type super-resolution methods—an extension of EXPRIT for non-Gaussian signals—can be considered as liner prediction of signals beyond both real and virtual antenna arrays.

Section 5 shows that the Capone spatial spectrum results from the harmonic averaging of maximum entropy estimates. Similarly, the MUSIC spectrum results from the harmonic averaging of Pisarenko estimates

Section 6 introduces a three-dimensional common space for angular super-resolution methods based on a three-dimensional space: the entropy of the correlation function beyond the antenna array aperture, the presence of harmonic averaging, and the use of multichannel signal prediction by subarrays. This space accommodates all super-resolution methods discussed in Section 3 to 5, along with any intermediate estimators. Two currently unoccupied boundary points in this space correspond to the multichannel variants of Capon and MUSIC. To illustrate the utility of this framework, we synthesize an estimator for one of these missing points—multichannel MUSIC. Additionally, we introduce a corresponding common space for non-Gaussian signal super-resolution.

In Section 7, we propose universal algorithms for estimating spatial super-resolution spectra via J orthogonalization. Grounded in the relationships between methods within a common space, this unified framework delivers simultaneous spectral estimates for various super-resolution techniques and any intermediate positions between them.

2. Linear Prediction

2.1. Maximum Entropy Linear Prediction

Linear prediction is used to extrapolate signals beyond the physical antenna aperture [8,9]. No specific behavior may be assumed for the correlation function outside the observation window, other than it exists and remains non-zero. In this case, the function can be extended to maximize randomness while ensuring the resulting function aligns with the measurements within the observation window. This approach corresponds to the Maximum Entropy method [2].

Consider a uniform linear array (ULA) with N antenna elements that are equally spaced with a distance is equal to a half wavelength. Assume that there are d far field narrowband signal sources located at ${\theta }_1,{\theta }_2,\dots ,{\theta }_d$. The received signals and noise are uncorrelated. Under these conditions, the output signal vector of the ULA with the number of snapshots t is as follows:

$$\mathbf{X}\left(t\right)={\left[x_1\left(t\right),x_2\left(t\right),\dots ,x_N\left(t\right)\right]}^T={\displaystyle \sum _{k=1}^d}\mathit{a}\left({\theta }_k\right)s_k\left(t\right)=\mathbf{A}\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right),$$

(1)

where $x_i$ is the output of the i’th array element; ${\theta }_k$ is the DOA of the k-th incident signal; $\mathbf{s}\left(t\right)$ is the signal vector; $\mathbf{n}\left(t\right)$ is the array additive noise vector; and $\mathbf{A}$ is the Vandermonde matrix, i’th column of which is the steering vector of the k’th signal:

$$\mathit{a}\left({\theta }_k\right)={\left[1,\mathrm{e}\mathrm{x}\mathrm{p}\left\{j2\pi {\displaystyle \frac{\mathsf{\Delta }}{\lambda }}\sin \left({\theta }_k\right),\dots ,j2\pi \left(N-1\right){\displaystyle \frac{\mathsf{\Delta }}{\lambda }}\sin \left({\theta }_k\right)\right\}\right]}^T,$$

(2)

where $\lambda$ is the wavelength of the incident signals, $\mathsf{\Delta }$ is the element spacing and is no more than half the wavelength, and the covariance matrix is as follows

$$\mathbf{R}=E\left\{\mathbf{X}\left(t\right){\mathbf{X}\left(t\right)}^{\mathit{H}}\right\}=\mathbf{A}\mathbf{S}{\mathbf{A}}^{\mathit{H}}+{\sigma }^2\mathbf{I}.$$

where $\mathbf{S}=E\left\{\mathbf{s}\left(t\right){\mathbf{s}\left(t\right)}^H\right\}$ is the signal covariance matrix, ${\sigma }^2$ is an unknown Gaussian noise variance.

The entropy of a stationary random process is given by the expression [10]

$$H={\displaystyle \underset{-\pi }{\overset{\pi }{\int }}}\mathrm{l}\mathrm{n}\left[P\left(\omega \right)\right]d\omega ,$$

(3)

where $\omega =2\pi {\displaystyle \frac{\mathsf{\Delta }}{\lambda }}\mathrm{s}\mathrm{i}\mathrm{n}\theta$ is the spatial frequency, and P(ω) is the spatial spectrum,

$$P\left(\omega \right)={\displaystyle \sum _{n=-\infty }^{n=\infty }}{r}_n{e}^{-j\omega n},$$

(4)

where r_n is a spatial correlation function.

Entropy can be parameterized directly through the determinant of the spatial correlation matrix [10],

$$H=\underset{N\to \infty }{\lim }{\displaystyle \frac{1}{2}}ln\left[\mathrm{d}\mathrm{e}\mathrm{t}{\mathbf{R}}_N\right],$$

(5)

where ${\mathbf{R}}_{\mathrm{N}}$ is a covariance matrix of N-element antenna array.

The maximum entropy estimate of the spatial spectrum will be

$$P_{ME}\left(\theta \right)={\displaystyle \frac{1}{{\mathbf{a}\left(\theta \right)}^H\mathbf{f}{\mathbf{f}}^H\mathbf{a}\left(\theta \right)}}={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H\mathbf{f}\right|}^2}}$$

(6)

where a(θ) = [1, e^jω, e^j2ω, …, e^j(N-1)ω]^T is the steering vector, f is the first column of ${\mathbf{R}}^{-1}$.

DOA estimates correspond to the peak locations of the spatial spectrum (6). The infinite correlation function in (5) with maximum entropy is uniquely determined by the vector f. Because this is an autoregressive process of order N, the elements of vector f serve as the corresponding autoregressive coefficients. Each subsequent value of the correlation function outside the observed window can be determined by the condition that the segment of the correlation function is orthogonal to the vector f. Consequently, estimating the spectrum using the autoregressive coefficients in (6) is equivalent to extrapolating the correlation function of the received signals beyond the antenna aperture using maximum entropy. In other words, maximum entropy spatial correlation function should be as close as possible to spatially white noise while matching the known points within the aperture.

2.2. Maximum to Minimum Entropy Linear Prediction

Maximum entropy approach assumes minimum knowledge about the behavior of the correlation function beyond the antenna array. However, this does not always fit radar surveillance, where some prior information about the correlation function beyond the observed segment is available. For far-field sources, wavefronts at the array can be treated as plane waves and extrapolated beyond the aperture. This approach constitutes a linear prediction framework; however, it optimizes for minimum entropy rather than maximum entropy.

The entropy of the spatial correlation function for Gaussian signals, as determined by expression (5), is uniquely related to the determinant of its covariance matrix. Consequently, the problem of minimizing this entropy is equivalent to minimizing the determinant of the covariance matrix for a stationary process extended to infinity beyond the observed region. In this case, due to the equidistance and stationarity of the spatial and temporal sample, the covariance matrix of signals is Toeplitz.

Extending the spatial covariance function$\left[r_0,r_1,\dots ,r_{N-1}\right]$ by one step corresponds to finding one new element $r_N$ of the Hermitian Toeplitz covariance matrix where we denote it as $\zeta$:

$${\mathbf{R}}_N=\left[\begin{array}{ccccc}{r}_0& r_1^{*}& \dots & r_{N-1}^{*}& {\zeta }^{*}\\ r_1& r_0& \dots & r_{N-2}^{*}& r_{N-1}^{*}\\ \dots & \dots & \dots & \dots & \dots \\ r_{N-1}& r_{N-2}& \dots & r_0& r_1^{*}\\ \zeta & r_{N-1}& \dots & r_1& r_0\end{array}\right],$$

(7)

ensuring the covariance matrix is semi-definite, $\det \left({\mathbf{R}}_N\right)\ge 0$.

By applying Sylvester's determinant identity [11], we can express the determinant of the Toeplitz covariance matrix ${\mathbf{R}}_N$ as a function of the predicted value:

$$\det \left({\mathbf{R}}_N\right)=D_N\left(\zeta \right)={\displaystyle \frac{D_{N-1}^2-{\left|D_{N-1,\left(1\right)}\left(\zeta \right)\right|}^2}{D_{N-2}}},$$

(8)

where

$$D_{N-1,\left(1\right)}\left(\zeta \right)=\det \left[\begin{array}{cccc}{r}_1& r_0& \dots & r_{N-2}^{*}\\ \dots & \dots & \dots & \dots \\ r_{N-1}& r_{N-2}& \dots & r_0\\ \zeta & r_{N-1}& \dots & r_1\end{array}\right].$$

As $\det \left({\mathbf{R}}_N\right)$ is a quadratic function of the predicted element $\zeta$, there are two values of $\zeta$, where the matrix determinant is equal to zero. These values define the boundaries for the solutions, ranging from minimum to maximum entropy as shown in Figure 1.

Since $D_{N-1}$ and $D_{N-2}$ are independent of $\zeta$, maximizing the determinant of the covariance matrix simplifies to minimizing the minor $D_{N-1,\left(1\right)}\left(\zeta \right)$. By applying the Sylvester identity to the determinant $D_{N-1,\left(1\right)}\left(\zeta \right)$, we transform (8) to

$$D_N\left(\zeta \right)={\displaystyle \frac{D_{N-1}^2-{\left|{\zeta D_{N-2}+D}_{N-1,\left(1\right)}\left(0\right)\right|}^2}{D_{N-2}}}.$$

(9)

This means, applying the maximum entropy condition yields the following equation

$${\zeta D_{N-2}+D}_{N-1,\left(1\right)}\left(0\right)=0,$$

(10)

whereas the minimum entropy condition gives the equation

$$D_{N-1}^2-{\left|{\zeta D_{N-2}+D}_{N-1,\left(1\right)}\left(0\right)\right|}^2=0.$$

(11)

Assuming D_N-2 > 0 and D_N (ζ) ≥ 0, we find a single solution of the maximum entropy and two solutions of the minimum entropy (D_N (ζ) = 0):

$${\zeta }_{\mathrm{m}\mathrm{a}\mathrm{x}}={\displaystyle \frac{D_{N-1,\left(1\right)}\left(0\right)}{D_{N-2}}},$$

(12)

$${\zeta }_{\mathrm{m}\mathrm{i}\mathrm{n}1}={\displaystyle \frac{D_{N-1,\left(1\right)}\left(0\right)-D_{N-1}}{D_{N-2}}},$$

(13)

$${\zeta }_{\mathrm{m}\mathrm{i}\mathrm{n}2}={\displaystyle \frac{D_{N-1,\left(1\right)}\left(0\right)+D_{N-1}}{D_{N-2}}}.$$

(14)

Minimum entropy roots may also be represented as

$${\zeta }_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{1,2}}={\zeta }_{\mathrm{m}\mathrm{a}\mathrm{x}}\pm {\displaystyle \frac{D_{N-1}}{D_{N-2}}}.$$

(15)

Expressions (12) and (15) determine the next predicted coefficient of the correlation function according to the criterion of maximum and minimum entropy, respectively. Because the extrapolated correlation function depends linearly on the first N known coefficients, the coefficients of a linear filter that generates the entire function can be directly determined. These filter coefficients define a spatial spectrum identical to that of the underlying infinite correlation function. Consequently, the infinite correlation function can be parameterized directly as a vector of linear prediction coefficients.

To derive the vector of linear prediction coefficients, we assume the maximum and minimum entropy correlation coefficients are given by (12) and (13), respectively. Accordingly, we partition the covariance matrix ${\mathbf{R}}_N$ into the following block form:

$${\mathbf{R}}_N=\left[\begin{array}{cc}{r}_0& {\mathbf{r}}_N^H\left(\zeta \right)\\ {\mathbf{r}}_{\mathbf{N}}\left(\zeta \right)& {\mathbf{R}}_{N-1}\end{array}\right].$$

(16)

Regardless of the variant of the criterion for extending covariance matrix, maximum or minimum entropy, the vector of coefficients of the linear prediction filter f is determined by the equation

$$\left[\begin{array}{cc}{r}_0& {\mathbf{r}}_N^H\left(\zeta \right)\\ {\mathbf{r}}_{\mathbf{N}}\left(\zeta \right)& {\mathbf{R}}_{N-1}\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{f}\end{array}\right]=\left[\begin{array}{c}{\rho }_N\\ 0\end{array}\right],$$

(17)

where ρ_N is a prediction error variance, which, for the minimum entropy solution, is assumed to be zero.

The columns of the covariance matrix ${\mathbf{R}}_{\mathrm{N}-1}$ in (17) are linearly independent and span the complete N-dimensional space containing the vector r_N(ζ). Therefore, regardless of the value of ζ, the vector r_N (ζ) can be expressed through the columns of the matrix${\mathbf{R}}_{\mathrm{N}-1}$,

$$\left[\begin{array}{cc}{\mathbf{r}}_{\mathbf{N}}\left(\zeta \right)& {\mathbf{R}}_{N-1}\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{f}\end{array}\right]=0,$$

(18)

hence

$${\mathbf{f}=-{\mathbf{R}}_{N-1}^{-1}\mathbf{r}}_{\mathbf{N}}\left(\zeta \right).$$

(19)

Substituting ζ $\in$ [ζ_min1, ζ_max] (or ζ $\in$ [ζ_max, ζ_min2]) defines linear prediction coefficient vectors spanning the maximum, minimum, and all intermediate entropy solutions.

For the maximum entropy solution $D_{N-1,\left(1\right)}\left({\zeta }_{max}\right)=0$. It means there is a vector f such that

$$\left[\begin{array}{cc}{\mathbf{r}}_{N-1}& {\mathbf{R}}_{N-2}\\ {\zeta }_{\mathrm{m}\mathrm{a}\mathrm{x}}& {\mathbf{r}}_{N-1}^H\overline{\mathbf{J}}\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{f}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right],$$

(20)

where $\overline{\mathbf{J}}$ is

$$\overline{\mathbf{J}}=\left[\begin{array}{ccc}& & 1\\ & ⋰& \\ 1& & \end{array}\right].$$

From (20), the coefficient vector for the maximum entropy prediction follows as:

$${\mathbf{f}=\mathbf{R}}_{N-2}^{-1}{\mathbf{r}}_{N-1},$$

(21)

and the next predicted correlation coefficient is

$${{\zeta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathbf{r}}_{N-1}^H\overline{\mathbf{J}}{\mathbf{R}}_{N-2}^{-1}{\mathbf{r}}_{N-1}.$$

(22)

Thus, the maximum entropy solution can be expressed in two ways: either by substituting the value of $\zeta$ from (12) into (19), or directly via expression (21). The resulting spatial spectra are identical because both approaches define the same extension of the covariance matrix.

Let us examine the minimum entropy solution and explore alternative methods for obtaining the minimum entropy linear prediction vector beyond expression (19).

In contrast to the strictly positive definite maximum entropy extension, the minimum entropy extension R_N of the covariance matrix R_N-1 is singular. At the first step of extending the covariance matrix, there are two variants of singular extension, which are defined in (13) and (14). Both options yield the same result since they generate equivalent spaces of the covariance matrix. Further singular extension of the matrix R_N is unique, as guaranteed by the infinite singular extension theorem [11]. In turn, the vector of linear prediction coefficients (19) and any variant of ζ (13) or (14) uniquely determine the covariance matrix extension and the corresponding minimum entropy spectrum.

The vector of the linear prediction coefficients is orthogonal to the space spanned by the rows of the covariance matrix R_N,

$$\left[\begin{array}{cc}{r}_0& {\mathbf{r}}_N^H\left({\zeta }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\\ {\mathbf{r}}_{\mathbf{N}}\left({\zeta }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)& {\mathbf{R}}_{N-1}\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{f}\end{array}\right]=\left[\begin{array}{c}0\\ \mathbf{0}\end{array}\right].$$

(23)

Consequently, the minimum-entropy linear prediction vector is the eigenvector corresponding to the zero eigenvalue of the extended covariance matrix R_N (ζ_min1,2).

Under the standard signal model, the covariance matrix R_{N -1} is full rank, whereas its minimum-entropy extension becomes singular. However, the minimum entropy criterion essentially assumes extrapolating signal wavefronts beyond the physical antenna aperture and does not require the original covariance matrix to be non-singular.

We consider a singular covariance matrix obtained by removing the uncorrelated noise component. The rank of the resulting matrix

$${\mathbf{C}}_{N-1}={\mathbf{R}}_{N-1}-{\sigma }^2\mathbf{I}=\mathbf{A}\mathbf{S}{\mathbf{A}}^H$$

(24)

equals the number of signal sources, $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathbf{A}\mathbf{S}{\mathbf{A}}^H\right)=d$.

Let us derive a recurrence formula for the singular matrix extension that yields the minimum-entropy linear prediction coefficients. Assume the noise-free covariance matrix have rank d<N. It follows that any d+1 columns of ${\mathbf{C}}_{N-1}$ are linearly dependent. Thus, for the first d+1 columns, there exist coefficients u₀, u₁, …, u_d, such that:

$$\left[\begin{array}{cccc}{c}_0& c_1^{*}& \dots & c_d^{*}\\ c_1& c_0& \dots & c_{d-1}^{*}\\ \dots & \dots & \dots & \dots \\ c_{N-1}& c_{N-2}& \dots & c_{N-d-1}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ \dots \\ u_{d+1}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \dots \\ 0\end{array}\right].$$

(25)

These coefficients are proportional to the cofactors of the elements in the first row of determinant D_d for matrix ${\mathbf{C}}_{N-1}$ [12]. Since appending a row preserves the matrix's rank, it follows that the first d+1 columns remain linearly dependent,

$$\left[\begin{array}{cccc}{c}_0& c_1^{*}& \dots & c_d^{*}\\ c_1& c_0& \dots & c_{d-1}^{*}\\ \dots & \dots & \dots & \dots \\ c_{N-1}& c_{N-2}& \dots & c_{N-d-1}\\ c_N& c_{N-1}& \dots & c_{N-d}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ \dots \\ u_{d+1}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \dots \\ 0\\ 0\end{array}\right],$$

(26)

where the singular extension element, ζ, is defined as

$$c_N=-{\displaystyle \frac{u_2}{u_1}}{c}_{N-1}-{\displaystyle \frac{u_3}{u_1}}{c}_{N-2}-{\displaystyle \frac{u_3}{u_1}}{c}_{N-2}-\dots -{\displaystyle \frac{u_{d+1}}{u_1}}{c}_{N-d}.$$

(27)

This defines a recurrence relation that recursively updates the spatial correlation function and its corresponding Toeplitz covariance matrix, preserving their inherent rank. The coefficients $u_i$ are the cofactors of the first-row elements of matrix C, though they can also be defined alternatively. Expression (25) implies that the forecast vector u is the right singular vector of the matrix formed by the first d+1 columns of С_N-1, associated with the zero singular value. By complementing u with zeros to match the larger dimension, the resulting vector lies in the null space of С_N-1, confirming it as an eigenvector of its zero eigenvalue:

$${\mathbf{C}}_{N-1}\left[\begin{array}{c}{\mathbf{u}}_{d+1}\\ 0_{N-d-1}\end{array}\right]=0.$$

(28)

Furthermore, because adding σ²I to С_N- simply shifts all its eigenvalues by σ² without altering its eigenvectors, we can express this as:

$${\mathbf{R}}_{N-1}\left[\begin{array}{c}{\mathbf{u}}_{d+1}\\ 0_{N-d-1}\end{array}\right]{=\sigma }^2\left[\begin{array}{c}{\mathbf{u}}_{d+1}\\ 0_{N-d-1}\end{array}\right].$$

(29)

Consequently, the zero-padded vector u acts as an eigenvector of the covariance matrix R_N-1, with an associated eigenvalue equal to the noise variance σ².

The entire correlation sequence can be generated using d initial values ($r_0,r_1,\dots ,r_d$) and the noise eigenvector, u, of the covariance matrix. While the Fourier transform of this infinite sequence yields the minimum entropy spatial spectrum, calculating it is unnecessary. Because the correlation sequence is recursive, its spectral characteristics are directly determined by the transfer function of a recursive linear filter represented by a polynomial:

$$H\left(z\right)={\displaystyle \frac{1}{u_1+u_2{z}^{-1}+u_3{z}^{-2}+\dots +u_{d+1}{z}^{-d}}}={\displaystyle \frac{1}{B\left(z\right)}},$$

(30)

where $z=\mathrm{e}\mathrm{x}\mathrm{p}\left(j2\pi {\displaystyle \frac{\mathsf{\Delta }}{\lambda }}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\right)$.

The spatial spectrum of signals is determined by the expression

$$P\left(\theta \right)={\displaystyle \frac{1}{{\left|{\displaystyle {\sum }_{k=1}^{d+1}}{u}_k\mathrm{e}\mathrm{x}\mathrm{p}\left(j2\pi \left(k-1\right)\frac{\mathsf{\Delta }}{\lambda }\sin \left(\theta \right)\right)\right|}^2}}={\displaystyle \frac{1}{{\left|{\mathbf{a}}^H\left(\theta \right)\mathbf{u}\right|}^2}}.$$

(31)

The spectral estimate in (31) is derived from the minimum entropy criterion of the extrapolated correlation function. Physically, this equates to extending the plane wavefronts of the received signals beyond the boundary of the antenna aperture. Because vector u in (31) is an eigenvector corresponding to the noise eigenvalue of the covariance matrix, this estimate also represents the Pisarenko spatial spectrum [4,13]. Therefore, the Pisarenko method extrapolates the correlation function beyond the observed region, much like the maximum entropy method. Unlike the latter, it relies on minimizing the entropy of the process's power spectral density. This corresponds to extrapolating the plane wavefronts of the sources beyond the antenna aperture.

3. Multi-Channel Linear Prediction

We define multi-channel linear prediction as a process where each prediction step extrapolates an antenna subarray rather than a single element. We demonstrate that the ESPRIT algorithm is fundamentally equivalent to this multi-channel linear prediction approach.

We begin with the classical formulation of the ESPRIT method, which we then reformulate within a linear prediction framework to interpret it as an extrapolation of source wavefronts beyond the antenna aperture. A distinctive feature of ESPRIT and related subspace methods [6,17,18,19,21,23] is that they require neither precise sensor calibration nor prior knowledge of the global array geometry. Instead, they rely on an array geometry composed of matched sensor pairs, where each pair shares an identical spatial displacement vector. The DOA estimates are then extracted directly from the generalized eigenvalues of a constructed matrix pencil.

In ESPRIT-like methods, it is convenient to consider the antenna array as consisting of two identical subarrays spatially displaced by the vector Δ. The antenna array receives narrowband signals with plane wavefronts arriving from external sources as presented in Figure 2. Let a_m = a(θ_m) denote the steering vector of the m-th source at the first subarray, corresponding to the angular coordinate θ_m. Because the subarrays are identical, the signals they receive from the same sources differ only by progressive phase shifts ϕ₁,..., ϕ_d. The output signals of the first and second subarrays are given by

$${\mathbf{X}}_1\left(t\right)=\left[\mathbf{a}\left({\theta }_1\right),\dots ,\mathbf{a}\left({\theta }_d\right)\right]\mathbf{s}\left(t\right)+{\mathbf{n}}_1\left(t\right)=\mathbf{A}\mathbf{s}\left(t\right)+{\mathbf{n}}_1\left(t\right),$$

(32)

$${\mathbf{X}}_2\left(t\right)=\left[\mathbf{a}\left({\theta }_1\right)e^{-j{\phi }_1},\dots ,\mathbf{a}\left({\theta }_d\right)\mathbf{a}{e}^{-j{\phi }_d}\right]\mathbf{s}\left(t\right)+{\mathbf{n}}_2\left(t\right)=\mathbf{A}\mathsf{\Phi }\mathbf{s}\left(t\right)+{\mathbf{n}}_2\left(t\right),$$

(33)

where n₁(t), n₂(t) are the additive noise vectors of the first and second subarrays; Φ = diag[exp{jϕ₁},..., exp{jϕ_d}], ϕ_m = 2|Δ| sinθ_m.

In generalized notation, the signals from the antenna array are expressed as follows:

$$\left(t\right)=\left[\begin{array}{c}{\mathbf{X}}_1\left(t\right)\\ {\mathbf{X}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}\mathbf{A}\\ \mathbf{A}\mathsf{\Phi }\end{array}\right]\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right)=\mathbf{G}\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right).$$

(34)

The source signals s(t) are zero-mean, temporary uncorrelated, complex Gaussian processes with a covariance matrix S = E[s(t)s^H(t)]. The noise n(t) also is modeled as a zero-mean, spatially and temporally white Gaussian process with variance σ². It is assumed to be uncorrelated with the source signals. The covariance matrix of the signals received by the antenna array is given by:

$$\mathbf{R}=E\left[\mathbf{X}\left(t\right){\mathbf{X}}^H\left(t\right)\right]=\mathbf{G}\mathbf{S}{\mathbf{G}}^H+{\sigma }^2.$$

(35)

We assume that matrices G and S have full rank rank(G) = rank(S) = d, implying that the source signals are not fully coherent.

In the ESPRIT method, the directions of arrival θ_m are determined from the non-zero generalized eigenvalues λ_m (m = 1, ..., d) of the matrix pencil {C₁₁, R₁₂}, which correspond to the non-zero solutions of the equation

$$\mathrm{d}\mathrm{e}\mathrm{t}\left[{\mathbf{C}}_{11}-\lambda {\mathbf{R}}_{12}\right]=0,$$

(36)

$${\theta }_m=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{{\lambda }_m}{2\pi \left|\mathsf{\Delta }\right|}}\right),$$

(37)

where ${\mathbf{C}}_{11}={\mathbf{R}}_{11}-{\sigma }^2$I is the covariance matrix of the first subarray without noise and R₁₂ is the cross-covariance matrix of the first and second subarrays:

$$\mathbf{R}=\left[\begin{array}{cc}{\mathbf{R}}_{11}& {\mathbf{R}}_{12}\\ {\mathbf{R}}_{21}& {\mathbf{R}}_{22}\end{array}\right],$$

(38)

$$\mathbf{C}=\mathbf{R}-{\sigma }^2{\mathbf{I}}_{2m}=\left[\begin{array}{cc}{\mathbf{C}}_{11}& {\mathbf{R}}_{12}\\ {\mathbf{R}}_{21}& {\mathbf{C}}_{22}\end{array}\right]=\left[\begin{array}{cc}\mathbf{A}\mathbf{S}{\mathbf{A}}^H& \mathbf{A}\mathbf{S}{{\mathsf{\Phi }}^H\mathbf{A}}^H\\ \mathbf{A}\mathsf{\Phi }\mathbf{S}{\mathbf{A}}^H& \mathbf{A}\mathsf{\Phi }\mathbf{S}{{\mathsf{\Phi }}^H\mathbf{A}}^H\end{array}\right]$$

(39)

We will now consider the ESPRIT method from the perspective of extrapolating the correlation function beyond the physical antenna array.

The antenna system described above — consisting of two identical subarrays spatially shifted by vector Δ — will be treated here as a modular array formed by two identical modules with shifted phase centers. First, we predict the cross-covariance matrix for the second and third (virtual) modules of the antenna array presented in Figure 3. Subsequent predictions are then generated recursively.

Extrapolating the wavefronts of signal sources in a modular antenna array is equivalent to enforcing the minimum entropy criterion. While a single step of uniform linear array (ULA) extrapolation predicts the signal or correlation function for a single virtual element, extrapolation in a modular array predicts either the signals across N spatial points or the cross-covariance matrix between the predicted and preceding subarrays.

Assuming spatial stationarity, the noise-free covariance matrix of the antenna array outputs is block-Toeplitz:

$$\mathbf{C}=\mathbf{R}-{\sigma }^2{\mathbf{I}}_{2m}=\left[\begin{array}{cc}{\mathbf{C}}_1& {\mathbf{C}}_2\\ {\mathbf{C}}_2^H& {\mathbf{C}}_1\end{array}\right].$$

(40)

Unlike conventional linear prediction which uses scalar-element vectors, multichannel linear prediction for wavefronts outside an antenna array utilizes matrix-element prediction vectors. This shift from scalars to matrices introduces additional degrees of freedom, which can be used to develop new methods for estimating matrix parameters and source angular coordinates under the minimum entropy criterion. We first introduce a block prediction matrix configuration to highlight the core principles of minimum-entropy multichannel prediction, and then demonstrate that ESPRIT estimates can be derived as a special case of this proposed model.

We now introduce the block matrix B,

$$\mathbf{B}=\left[\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{B}}_2\end{array}\right],$$

where B₁, B₂ are N×N matrices satisfy the condition

$$\mathbf{C}\mathbf{B}=\left[\begin{array}{cc}{\mathbf{C}}_1& {\mathbf{C}}_2\\ {\mathbf{C}}_2^H& {\mathbf{C}}_1\end{array}\right]\left[\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{B}}_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$

(41)

Note that rank(C₁) = rank(C₂) = d.

The block matrix B can be viewed as a generalized eigenvector of the matrix C corresponding to a zero eigenvalue. Because

$$\dim \left(\ker \left(\mathbf{C}\right)\right)=2N-d>N,$$

(42)

these vectors are not unique, which allows for some flexibility in its estimation.

On the other hand, equation (36) defines the generalized eigenvalues and eigenvectors of the matrix pencil {C₁₁, R₁₂} (denoted here as {C₁, C₂}). This can be expressed as

$$\left({\mathbf{C}}_1-\lambda {\mathbf{C}}_2\right)\mathbf{u}=\mathbf{0},$$

(43)

where λ and u are the generalized eigenvalue and corresponding generalized eigenvector, respectively. Equation (43) can be rewritten as follows:

$$\left[\begin{array}{cc}{\mathbf{C}}_1& {\mathbf{C}}_2\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ -\lambda \mathbf{u}\end{array}\right]=0.$$

(44)

A vector u satisfying (43) or (44) is a generalized eigenvector of the pencil {C₁, C₂}, meaning the vector [u^T, –λu^T]^T is both a right singular vector of [C₁, C₂] (zero singular value) and an eigenvector of C (zero eigenvalue), thus satisfying condition (41) for the columns of B. The block matrix B can thus be defined as

$${\mathrm{B}}_1=\mathrm{U},{\mathrm{B}}_2=\text{–}\mathrm{U}\mathsf{\Lambda },\mathbf{B}=\left[\begin{array}{c}\mathbf{U}\\ -\mathbf{U}\mathsf{\Lambda }\end{array}\right],$$

(45)

where U = [u₁, u₂, ..., u_N] and Λ = diag{λ_m}. Here, u_m and λ_m (for m=1, ..., N) denote the m-th generalized eigenvector and the corresponding generalized eigenvalue of the pencil {C₁, C₂}.

Order the eigenvalues λ_m so that λ_d+1 = λ_d+2 = … = λ_N = 0. Consequently, the vectors u_d+1, u_d+2, ..., u_N are also eigenvectors of matrix C₁ corresponding to its zero eigenvalues, while u₁, u₂, ..., u_d are proportional to the columns of the matrix $\mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}{\mathsf{\Phi }}^H$.

Therefore, matrix B can be also expressed as:

$$\mathbf{B}=\left[\begin{array}{cc}\mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}& {\mathbf{U}}_{d+1\dots N}\\ \mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}{\mathsf{\Phi }}^H& 0\end{array}\right],$$

(46)

where U_d+1...N is a matrix consisting of the eigenvectors associated with the zero eigenvalues of the subarray covariance matrix. Because matrix B is not uniquely defined, we present one possible definition to demonstrate the existence of a matrix that satisfies the above conditions. Practical methods for estimating B are discussed below.

The block matrix B defined in (45) acts as the multichannel linear predictor for the subarray, shifted by the vector Δ. Furthermore, the cross-covariance matrix C₃ between the second and third subarrays is defined by the following expression:

$${\mathbf{C}}_3\mathbf{U}\mathsf{\Lambda }={\mathbf{C}}_2\mathbf{U},$$

(47)

and for an arbitrary (k+1)-th step of predicting, respectively

$${\mathbf{C}}_{k+1}\mathbf{U}\mathsf{\Lambda }={\mathbf{C}}_k\mathbf{U}.$$

(48)

Equation (48) is a multichannel autoregressive model, extending the standard autoregressive formulation to account for multiple channels.

For this antenna array configuration, the steering vector corresponding to an incoming signal from direction θ is structured as follows:

$$\mathbf{a}\left(\theta \right)=\left[\begin{array}{c}{\mathbf{a}}_1\left(\theta \right)\\ \gamma {\mathbf{a}}_1\left(\theta \right)\end{array}\right],$$

(49)

where a₁(θ) is the amplitude-phase distribution vector of the reference signal received from direction θ at the aperture of the first subarray, and γ = exp(j2π|Δ| sinθ).

Following the scalar case, the multichannel spatial spectrum is determined by its block-matrix prediction filter:

$$P\left(\theta \right)={\displaystyle \frac{1}{{‖\mathbf{a}{\left(\theta \right)}^H\mathbf{B}‖}^2}}.$$

(50)

Without specifying the structure of matrix B, and relying solely on the fact that its columns span the noise subspace of the covariance matrix R (or the null space of matrix C), we can find that expression (50) represents the MUSIC spatial spectrum estimate. However, by exploiting the specific structure of the eigenvector matrix B, the extremum condition for functional (50) can be reformulated as

$${\mathbf{a}}^H\left(\theta \right)\mathbf{B}=0,$$

(51)

from where, considering (45) and (49), we obtain

$${\mathbf{a}}^H\left(\theta \right)\mathbf{U}=\gamma {\mathbf{a}}^H\left(\theta \right)\mathbf{U}\mathsf{\Lambda },$$

(52)

$${\gamma }_m={\lambda }_m^{*}=\exp \left(j2\pi \left|\mathsf{\Delta }\right|\sin \left({\theta }_m\right)\right),m=1,\dots ,d.$$

Thus, the ESPRIT spatial spectrum estimate can be interpreted as a multichannel autoregressive prediction of signals beyond the antenna aperture based on the minimum entropy criterion. Unlike the single-channel minimum entropy prediction in Pisarenko case, ESPRIT predicts the plane wavefronts of the sources across the aperture of the next N-element subarray.

Similar to the single-channel prediction, the multi-channel prediction can be optimized using not only the minimum entropy criterion, but also the maximum entropy criterion. Below we derive the maximum entropy spatial spectrum estimate for multichannel subarray prediction.

The signal and covariance matrix models for the antenna array remain unchanged, corresponding to (34) and (35). In the single-channel case, the maximum entropy solution reduces to a single row (or column) of the signal covariance matrix. For the multichannel case, this scalar covariance structure is extended, with the standard covariance matrix replaced by a corresponding block matrix:

$${\left[\begin{array}{cc}{\mathbf{R}}_{11}& {\mathbf{R}}_{12}\\ {\mathbf{R}}_{12}^H& {\mathbf{R}}_{11}\end{array}\right]}^{-1}=\left[\begin{array}{cc}{\mathbf{K}}^{-1}& {\mathbf{K}}^{-1}{\mathbf{R}}_{12}{\mathbf{R}}_{11}^{-1}\\ -{\mathbf{R}}_{11}^{-1}{\mathbf{R}}_{12}^H{\mathbf{K}}^{-1}& {\mathbf{R}}_{11}^{-1}+{\mathbf{R}}_{11}^{-1}{\mathbf{R}}_{12}^H{\mathbf{K}}^{-1}{\mathbf{R}}_{12}{\mathbf{R}}_{11}^{-1}\end{array}\right],$$

(53)

where

$$\mathbf{K}={{\mathbf{R}}_{11}-\mathbf{R}}_{11}^{-1}{\mathbf{R}}_{12}^H.$$

Matrix K normalizes the first column of the inverse block matrix (53) so that its first block becomes the identity matrix. Right-multiplying the first block column of (53) by K yields the block vector of the multichannel autoregressive maximum entropy prediction:

$$\mathbf{F}=\left[\begin{array}{c}\mathbf{I}\\ {-\mathbf{R}}_{11}^{-1}{\mathbf{R}}_{12}^H\end{array}\right].$$

(54)

The multichannel estimate of the maximum entropy spatial spectrum is given by:

$$P_{ESPME}\left(\theta \right)={\displaystyle \frac{1}{{‖\mathbf{a}{\left(\theta \right)}^H\mathbf{F}‖}^2}}={\displaystyle \frac{1}{{‖{\mathbf{a}}_1{\left(\theta \right)}^H-{\mathbf{a}}_2{\left(\theta \right)}^H{\mathbf{R}}_{11}^{-1}{\mathbf{R}}_{12}^H‖}^2}},$$

(55)

where a₁(θ) и a₂(θ) are the steering vectors for the first and second subarrays, respectively.

The angular coordinate estimates derived from the spatial spectrum (55) correspond to the zeros of the function

$${‖{\mathbf{a}}_1{\left(\theta \right)}^H-{\mathbf{a}}_2{\left(\theta \right)}^H{\mathbf{R}}_{11}^{-1}{\mathbf{R}}_{12}^H‖}^2=0.$$

(56)

Due to the specific configuration of the antenna array, the amplitude-phase distribution vectors are collinear, i.e., a₁(θ) = λa₂(θ), where λ = exp(-j2π|Δ| sinθ). Consequently, equation (56) simplifies to:

$${\mathbf{R}}_{12}{\mathbf{R}}_{11}^{-1}{\mathbf{a}}_1\left(\theta \right)=\lambda {\mathbf{a}}_1\left(\theta \right).$$

(57)

The equation provides exact solutions when the correlation function of the received signals follows a maximum-entropy autoregressive model. Specifically, these solutions are the eigenvalues of the matrix ${\mathbf{R}}_{12}{\mathbf{R}}_{11}^{-1}$, which define the angular coordinates of the signal sources.

Consequently, the approximate ESPRIT solutions are obtained by solving the standard eigenvalue problem for the matrix ${\mathbf{R}}_{12}{\mathbf{R}}_{11}^{-1}$. These solutions effectively predict the signal wavefronts beyond the physical antenna aperture while maximizing the entropy of the block correlation function.

4. Multi-Channel Linear Prediction of Non-Gaussian Signals

The Virtual-ESPRIT Algorithm (VESPA) [7,19,20] extends matrix pencil methods by utilizing fourth-order cumulants for non-Gaussian signals. Like standard ESPRIT, VESPA extracts source angles from the generalized eigenvalues of a matrix pencil. However, while standard approaches require a second physical antenna array to construct a cross-covariance matrix, VESPA replaces it with a virtual array. Cross-correlations are estimated directly from the cumulants of the single physical array, eliminating the need for additional subarray.

In this section, we reframe VESPA-type methods. We transition from their classical, cumulant-matrix definition to their interpretation as linear signal prediction beyond the physical antenna aperture.

Unlike the original ESPRIT method, VESPA uses a single N-element antenna array of arbitrary configuration presented in Figure 4. Two channels of this array must be identical and spatially shifted by a known vector Δ.

Calculating the angular coordinates of the signal sources using the matrix pencil method requires estimating both the spatial covariance matrix of the main antenna array and the cross-covariance matrix between the main array and its spatial shifted version by the vector Δ. This shifted copy is physically absent, serving instead as the virtual array illustrated in Figure 4.

Cumulant matrix of the physical antenna array for non-Gaussian sources takes the following form [7]:

$${\tilde{\mathbf{R}}}_{XX}=\mathbf{A}\mathsf{\Gamma }{\mathbf{A}}^H,$$

(58)

where A is the array steering matrix of the signal wavefronts on the real antenna aperture and Γ is the diagonal matrix of the source signals' fourth-order cumulants, $\mathsf{\Gamma }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{\gamma }_k^{\left(4\right)}\right\}$.

The mutual fourth-order cumulant matrix of the real and virtual antenna arrays is defined as follows:

$${\tilde{\mathbf{R}}}_{XV}={\displaystyle \sum _{k=1}^d}{\gamma }_k^{\left(4\right)}{e}^{j{\phi }_k}{\mathbf{a}}_k{\mathbf{a}}_k^H=\mathbf{A}\mathsf{\Gamma }{\mathsf{\Phi }}^H{\mathbf{A}}^H,$$

(59)

where ${\phi }_k=2\pi \left|\Delta \right|\sin {\theta }_k$ and ${\theta }_k$ is the direction of arrival of the k-th source signal.

By analyzing the structure of cumulant matrices (65) and (69), we show that the generalized eigenvalues of the matrix pencil $\left\{{\tilde{\mathbf{R}}}_{XX},{\tilde{\mathbf{R}}}_{XV}\right\}$ directly yield the values ${\xi }_k=e^{j{\phi }_k}$ corresponding to the signal bearings. This behavior mirrors the {C_XX, R_XV} pencil that the ESPRIT method would employ if a second antenna array existed. This can be verified through the pencil equation. Since

$${\tilde{\mathbf{R}}}_{XX}-\xi {\tilde{\mathbf{R}}}_{XV}=\mathbf{A}\mathsf{\Gamma }{\mathbf{A}}^H-\xi \mathbf{A}\mathsf{\Gamma }{\mathsf{\Phi }}^H{\mathbf{A}}^H=\mathbf{A}\left(\mathsf{\Gamma }{\mathbf{A}}^H-\xi \mathbf{A}\mathsf{\Gamma }{\mathsf{\Phi }}^H\right){\mathbf{A}}^H,$$

(60)

and both Γ and Ф are diagonal matrices, the rank of the matrix ${\tilde{\mathbf{R}}}_{XX}-\xi {\tilde{\mathbf{R}}}_{XV}$ jumps precisely when the parameters ξ_k equal the diagonal elements of Ф.

Thus, higher-order statistics enable the prediction of source wavefronts on a virtual aperture without constraining the antenna array configuration. Unlike the ESPRIT method, which requires N pairs of identical sensors, only a single identical pair is necessary here. Furthermore, the inherent noise-filtering properties of fourth-order cumulants eliminate the need to pre-process the covariance matrix to remove spatially uncorrelated noise. This suppression applies to both white and time-correlated noise. In practice, however, effective noise reduction depends on averaging a sufficient number of signal samples.

We employ higher-order statistics to predict the field distribution across the aperture of a second identical antenna array. This represents the first step in extending the plane wavefronts of signal sources beyond the existing antenna aperture. As a result, we obtain two cumulant matrices, ${\tilde{\mathbf{R}}}_{XX}$ and ${\tilde{\mathbf{R}}}_{XV}.$These matrices can be represented as sums of dyadic wavefront matrices originating from sources at two identical, spatially offset apertures. Our objective is to extend the source wavefronts further while maintaining a fixed rank for the cumulant matrices. Since matrices derived from fourth-order cumulants differ structurally from conventional covariance matrices, we analyze this process in detail.

The full cumulant matrix of the signals received by the antenna array, which comprises both real and virtual subarrays, takes the form:

$$\tilde{\mathbf{R}}=\left[\begin{array}{cc}{\tilde{\mathbf{R}}}_{XX}& {\tilde{\mathbf{R}}}_{XV}\\ {\mathbf{R}}_{XV}^H& {\tilde{\mathbf{R}}}_{XX}\end{array}\right].$$

(61)

The matrix has a block-Toeplitz structure and can be defined directly by the source wavefronts:

$$\tilde{\mathbf{R}}=\left[\begin{array}{cc}\mathbf{A}\mathsf{\Gamma }{\mathbf{A}}^H& \mathbf{A}\mathsf{\Gamma }{\mathsf{\Phi }}^H{\mathbf{A}}^H\\ \mathbf{A}\mathsf{\Phi }\mathsf{\Gamma }{\mathbf{A}}^H& \mathbf{A}\mathsf{\Gamma }{\mathbf{A}}^H\end{array}\right].$$

(62)

Note that the complete antenna array's cumulant matrix adopts form (62) only when all submatrices are estimated using fourth-order cumulants. If the matrix of the real antenna array is replaced by second-order statistics—which are available for the real array but not for the virtual one—then the generalized eigenvalues of the pencil $\left\{{\tilde{\mathbf{R}}}_{XX},{\tilde{\mathbf{R}}}_{XV}\right\}$ will no longer determine the bearings of the signal sources. Consequently, it becomes impossible to extrapolate the sources' plane wavefronts beyond the real and virtual arrays using this matrix combination.

Formula (5) is inapplicable because it relies on the assumption of Gaussian statistics. For non-Gaussian signals, maximizing or minimizing process entropy is no longer equivalent to optimizing the determinant of the covariance matrix. Consequently, to extend the signal beyond the antenna aperture, we return to the criterion for plane wavefront extension with a fixed number of non-Gaussian sources (the cumulant matrix rank). Formally, this requires finding a block-Toeplitz extension of the cumulant matrix (61) that retains the rank of the original matrix.

By analogy with the extension of covariance matrices, we introduce the block matrix

$$\mathbf{B}=\left[\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{B}}_2\end{array}\right]$$

whose blocks are square matrices B₁ and B₂ of dimension N×N, such that

$$\tilde{\mathbf{R}}\mathbf{B}=\left[\begin{array}{cc}{\tilde{\mathbf{R}}}_{XX}& {\tilde{\mathbf{R}}}_{XV}\\ {\mathbf{R}}_{XV}^H& {\tilde{\mathbf{R}}}_{XX}\end{array}\right]\left[\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{B}}_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$

(63)

Since the projected cumulant matrix ${\tilde{\mathbf{R}}}_{\infty }$ is block Toeplitz and its rank is fixed, ${\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\widetilde{(\mathbf{R}}}_{\infty })=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\tilde{\mathbf{R}}\right)={\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\tilde{\mathbf{R}}}_{XX})=d$, the block eigenvector of the cumulant matrix $\tilde{\mathbf{R}}$ uniquely determines all block eigenvectors spanning the null space of the extended cumulant matrix, which can be expressed as follows:

$$\left[\begin{array}{cccc}{\tilde{\mathbf{R}}}_{XX}& {\tilde{\mathbf{R}}}_{XV}& {\tilde{\mathbf{R}}}_{XZ}& \dots \\ {\mathbf{R}}_{XV}^H& {\tilde{\mathbf{R}}}_{XX}& {\tilde{\mathbf{R}}}_{XV}& \dots \\ {\mathbf{R}}_{XZ}^H& {\mathbf{R}}_{XV}^H& {\tilde{\mathbf{R}}}_{XX}& \dots \\ \dots & \dots & \dots & \dots \end{array}\right]\left[\begin{array}{cccc}{\mathbf{B}}_1& 0& 0& \dots \\ {\mathbf{B}}_2& {\mathbf{B}}_1& 0& \dots \\ 0& {\mathbf{B}}_2& {\mathbf{B}}_1& \dots \\ \dots & \dots & \dots & \dots \end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ \dots \end{array}\right],$$

(64)

where ${\tilde{\mathbf{R}}}_{XZ}$ denotes the subsequent predicted mutual cumulant matrix between the physical antenna array X and the predicted virtual array Z. The antenna array configuration is presented in Figure 5.

The validity of (64) can be demonstrated by isolating any 2N×2N diagonal block in ${\tilde{\mathbf{R}}}_{\infty }$, which corresponds to two adjacent virtual antenna arrays. Because ${\tilde{\mathbf{R}}}_{\infty }$ is block Toeplitz, this diagonal submatrix is identical to $\tilde{\mathbf{R}}$ and thus shares its block eigenvector, В. Since the rank of ${\tilde{\mathbf{R}}}_{\infty }$ does not exceed that of $\tilde{\mathbf{R}}$, any submatrix row above or below this block is a linear combination of $\tilde{\mathbf{R}}$. Consequently, В, padded with appropriately dimensioned zero blocks, forms a block eigenvector of ${\tilde{\mathbf{R}}}_{\infty }$ associated with its zero eigenvalue.

Expression (64) shows that the block vector B acts as a predictor for the cumulant matrices of the real and virtual subarrays. Specifically, the cumulant matrix ${\tilde{\mathbf{R}}}_{XZ}$ can be expressed as the solution to:

$${\tilde{\mathbf{R}}}_{XV}{\mathbf{B}}_1=-{\tilde{\mathbf{R}}}_{XZ}{\mathbf{B}}_2.$$

(65)

Alternatively, the matrix pencil $\left\{{\tilde{\mathbf{R}}}_{XX},{\tilde{\mathbf{R}}}_{XV}\right\}$ can be expressed in the form

$$\left[\begin{array}{cc}{\tilde{\mathbf{R}}}_{XX}& {\tilde{\mathbf{R}}}_{XV}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ -\xi \mathbf{u}\end{array}\right]=0,$$

(66)

where u is the generalized eigenvector of the pencil, and ξ is the generalized eigenvalue.

Comparing (63) and (66) shows the vectors ${\left[{\mathbf{u}}^T,-\xi {\mathbf{u}}^T\right]}^T$ in (66) lie in the column space of B. Following the logic of the matrix pencil method, we can determine the components of the block matrix B as follows:

$${\mathbf{B}}_1=\mathbf{U};{\mathbf{B}}_2=-\mathbf{U}\mathsf{\Xi },$$

(67)

where $\mathbf{U}=\left[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_N\right]$ and $\mathsf{\Xi }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}{\left\{{\xi }_k\right\}}_{k=1,\dots ,N}$.

To confirm the existence of a block vector B in the form (67), we provide an example for the case where the true cumulant matrix has the form (62). Of the N generalized eigenvalues of the pencil $\left\{{\tilde{\mathbf{R}}}_{XX},{\tilde{\mathbf{R}}}_{XV}\right\}$, exactly $\left(N-d\right)$values of ξ_k are zero, while ${\xi }_{d+1},\dots ,{\xi }_N$ are assumed to be nonzero. Consequently, the eigenvectors ${\mathbf{u}}_{d+1},{\mathbf{u}}_{d+2},\dots ,{\mathbf{u}}_N$ span the null space of ${\tilde{\mathbf{R}}}_{XX}.$ The remaining eigenvectors ${\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_d$ lie in the column space of A(A^HA)^-1Ф^H. This demonstrates that the block vector B can take the form (46), consistent with the ESPRIT method.

The predicted cross-cumulant matrix of the X and Z subarrays is defined by the following equation:

$${\tilde{\mathbf{R}}}_{XZ}\mathsf{\Xi }={\tilde{\mathbf{R}}}_{XV}.$$

(68)

At an arbitrary step (k+1), the cumulant matrix for the (k+1)-th antenna array satisfies the equation:

$${\tilde{\mathbf{R}}}_{k+1}\mathbf{U}\mathsf{\Xi }={\tilde{\mathbf{R}}}_k\mathbf{U}.$$

(69)

Equation (68) defines a zero-residual multichannel autoregressive model with matrix coefficients ${\mathbf{B}}_1=\mathbf{U}$ and ${\mathbf{B}}_2=-\mathbf{U}\mathsf{\Xi }$. Consequently, the VESPA virtual aperture expansion method can be interpreted as an extrapolation of plane wavefronts beyond both the real and the virtual subarrays. This extrapolation is driven by a block multichannel, minimum-entropy autoregressive process.

5. Averaging Linear Prediction Estimates

Previous sections established that the Pisarenko, ESPRIT, and VESPA algorithms can be interpreted as techniques for extrapolating the signal correlation function or wavefronts beyond the antenna aperture. Yet, certain angular super-resolution methods—such as Capon [3] and MUSIC [5] (and its variants [18,22,23,24,25,26])—do not inherently extend a single correlation function beyond the observed domain. We will demonstrate that these techniques can still be considered as liner prediction methods, provided the corresponding spatial spectrum estimates are averaged over subarrays.

The following expressions define the spatial spectrum estimates for the Capon and maximum entropy methods, respectively:

$$P_{CP}\left(\theta \right)={\displaystyle \frac{1}{{\mathbf{a}\left(\theta \right)}^H{\mathbf{R}}^{-1}\mathbf{a}\left(\theta \right)}},$$

(70)

$$P_{ME}\left(\theta \right)={\displaystyle \frac{1_1^T{\mathbf{R}}^{-1}{1}_1}{{\left|{\mathbf{a}\left(\theta \right)}^H\mathbf{f}\right|}^2}},$$

(71)

where R is the covariance matrix of the received signals, f is the vector of autoregressive parameters satisfying the maximum entropy criterion, a(θ) is the steering vector, and 1₁ = [1, 0, ..., 0]^T.

The vector of autoregressive parameters, f, represents the first column of the inverse covariance matrix ${\mathrm{R}}^{-1}$. Therefore, expression (71) can be rewritten as:

$$P_{ME}\left(\theta \right)={\displaystyle \frac{1_1^T{\mathbf{R}}^{-1}{1}_1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{R}}^{-1}{1}_1\right|}^2}}.$$

(72)

To demonstrate the relationship between the Capon method and maximum entropy autoregressive prediction, we use a triangular factorization of the inverse covariance matrix,

$${\mathrm{R}}^{-1}=\mathrm{L}{\mathrm{L}}^H,$$

(73)

where L is a lower triangular matrix.

The Capon and maximum entropy spectral estimates are then given by the following expressions:

$$P_{CP}\left(\theta \right)={\displaystyle \frac{1}{{\mathbf{a}\left(\theta \right)}^H\mathrm{L}{\mathrm{L}}^H\mathbf{a}\left(\theta \right)}}={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H\mathrm{L}\right|}^2}}={\displaystyle \frac{1}{{\displaystyle \sum _{k=1}^N}{\left|{\mathbf{a}\left(\theta \right)}^H{l}_{*k}\right|}^2}},$$

(74)

$$P_{ME}\left(\theta \right)={\displaystyle \frac{1_1^T\mathrm{L}{\mathrm{L}}^H{1}_1}{{\left|{\mathbf{a}\left(\theta \right)}^H\mathrm{L}{\mathrm{L}}^H{1}_1\right|}^2}}={\displaystyle \frac{{\left|l_{11}\right|}^2}{{\left|{\mathbf{a}\left(\theta \right)}^H{l}_{*1}{l}_{11}\right|}^2}}={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{l}_{*1}\right|}^2}},$$

(75)

where l_*k denotes the k-th column of the matrix L.

A comparison of (74) and (75) reveals the relationship between the Capon and maximum entropy estimates:

$$P_{CP}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^N}\frac{1}{P_{ME,k}\left(\theta \right)}}},$$

(76)

where P_ME,k(θ) is the autoregressive estimate of the maximum entropy spectrum of order k.

The estimate P_ME,k(θ) corresponds to the standard maximum entropy spectrum of a subarray formed by the first k elements of the original N-element antenna:

$$P_{ME,k}\left(\theta \right)={\displaystyle \frac{{\mathbf{1}}_1^T{\mathbf{R}}_k^{-1}{\mathbf{1}}_1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{R}}_k^{-1}{\mathbf{1}}_1\right|}^2}},$$

(77)

where ${\mathbf{R}}_k^{-1}$ is the inverse covariance matrix of the first k elements subarray.

Therefore, the Capon spectrum estimate represents the harmonic mean of the maximum entropy spectral estimates across all model orders from 1 to N. This explains the lower resolution of the Capon method. While the denominator in the maximum entropy spatial spectrum (75) corresponds to the adaptive pattern of a full N-element antenna array, the denominator in the Capon spatial spectrum (74) represents the sum of adaptive patterns across arrays of sizes 1, 2, …, N. The lower-order patterns smooth the aggregate adaptive pattern, broadening the main beam of the Capon spectrum.

Figure 6 compares the maximum entropy and Capon estimates, interpreting both as methods for extrapolating the correlation function beyond the observed window. Here, F_∞ denotes the Fourier transform of an infinite sequence. A k-order recursive filter represents a linear feedback filter that generates a k-order autoregressive process. The filter coefficients are defined as the elements in the first column of the inverse covariance matrix of the corresponding order, while the initialization values are the corresponding correlation coefficients.

Generating an infinite correlation sequence and computing its spatial spectrum via the Fourier transform is mathematically equivalent to convolving the autoregressive parameters with steering vectors across all directions in according to (74) and (75). Figure 6 illustrates the conceptual interpretation of the maximum entropy and Capon methods regarding correlation function extrapolation, rather than their specific computational steps.

Figure 7 illustrates the relationship between the maximum entropy and Capon spectra defined in (76). The Capon spectrum in Figure 7(i) is obtained by averaging the maximum entropy spectra across antenna subarrays with 1 to 16 elements.

The graphs in Figure 7 illustrate the relationship between the maximum entropy and Capon spectra, as defined in (76), showing that the Capon spectrum presented in Figure 7(i) is derived by averaging maximum entropy spectra across antenna subarrays with elements ranging from 1 to 16. While this averaging process reduces sensitivity to antenna array calibration errors, it degrades spatial resolution, failing to resolve sources at coordinates 0.15 and 0.22.

Next, we consider MUSIC method, a direct generalization of the Pisarenko approach. Their respective spatial spectrum estimates are defined as:

$$P_{MUSIC}\left(\theta \right)={\displaystyle \frac{1}{{\mathbf{a}\left(\theta \right)}^H{\mathrm{E}}_{N-d}{\mathrm{E}}_{N-d}^H\mathbf{a}\left(\theta \right)}},$$

(78)

$$P_{PSR}\left(\theta \right)={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{e}}_1\right|}^2}},$$

(79)

where E_N-d = [e₁, e₂, …, e_N-d] is an $N\times (N-d)$ matrix whose columns are the noise eigenvectors corresponding to the repeated noise eigenvalue σ² of the covariance matrix R.

The MUSIC spectrum can be also expressed as

$$P_{MUSIC}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^{N-d}}{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{e}}_k\right|}^2}}.$$

(80)

While both MUSIC and Pisarenko are subspace methods, viewing MUSIC as an extrapolation of the autocorrelation function requires distinct analyses for single repeated noise eigenvalues and unique noise eigenvalues.

First, we consider the case where the covariance matrix has a Toeplitz structure, and its spectrum contains a single repeated noise eigenvalue. We show that, under these conditions, the Pisarenko and MUSIC methods yield equivalent estimates.

The noise subspace projection matrix ${\mathbf{E}}_{N-d}{\mathbf{E}}_{N-d}^H$ of the covariance matrix in (85) can be represented as:

$${\mathrm{E}}_{N-d}{\mathrm{E}}_{N-d}^H={\mathrm{E}}_{N-d}{{\mathrm{Q}}^H\mathrm{Q}\mathrm{E}}_{N-d}^H=\mathbf{W}{\mathbf{W}}^H,$$

(81)

where Q be an (N-d)×(N-d) unitary matrix (Q^HQ = I) that transforms the noise eigenvectors of the covariance matrix into a lower triangular matrix W presented in Figure 8. Although the columns of W are not mutually orthogonal, they still span the noise subspace of the covariance matrix.

Since a unitary transformation of the eigenvectors preserves the underlying subspace, the MUSIC spectrum can be expressed as:

$$P_{MUSIC}\left(\theta \right)={\displaystyle \frac{1}{{\mathbf{a}\left(\theta \right)}^H\mathbf{W}{\mathbf{W}}^H\mathbf{a}\left(\theta \right)}}={\displaystyle \frac{1}{{\displaystyle \sum _{k=1}^{N-d}}{\left|{\mathbf{a}\left(\theta \right)}^H{w}_k\right|}^2}},$$

(82)

where w_k is the k-th column of the matrix W.

Since W is triangular, the subvector of non-zero elements in column w_k, (omitting the first) serves as an eigenvector for the subarray of the last (N-k+1) antenna elements. This relationship is written as:

$${\mathbf{R}}^{\left(k\right)}{\mathbf{W}}^{\left(k\right)}={\sigma }^2{\mathbf{W}}^{\left(k\right)},$$

(83)

where R^(k) and W^(k) are the lower-right submatrices of R and W, with dimensions (N-k+1)×(N-k+1) and (N-k+1)×(N-d-k+1), respectively,

$${\mathbf{W}}^{\left(k\right)}=\left[\begin{array}{cc}{0}_{k-1}& 0\\ 0& {\mathbf{I}}_{N-k+1}\end{array}\right]\mathbf{W}\left[\begin{array}{cc}{0}_{k-1}& 0\\ 0& {\mathbf{I}}_{N-d-k+1}\end{array}\right],$$

(84)

$${\mathbf{R}}^{\left(k\right)}=\left[\begin{array}{cc}{0}_{k-1}& 0\\ 0& {\mathbf{I}}_{N-k+1}\end{array}\right]\mathbf{R}\left[\begin{array}{cc}{0}_{k-1}& 0\\ 0& {\mathbf{I}}_{N-k+1}\end{array}\right].$$

(85)

As established above, the noise eigenvector of the Toeplitz covariance matrix defines an autoregressive extrapolation of the correlation function beyond the antenna aperture. Consequently, each noise eigenvector w_k of dimension (N-k+1) defines a minimum-entropy extrapolation for the (N-k+1)-element subarray. The autoregressive equation describing this extension is given by:

$$r_{s+1}^{\left(k\right)}=-{\displaystyle \frac{1}{w_{N-k+1}^{\left(k\right)}}}{\displaystyle \sum _{m=1}^{N-k}}{w}_m^{\left(k\right)}{r}_{s+i-N+k}^{\left(k\right)},$$

(86)

where $r_{s+1}^{\left(k\right)}$is the (s+1)-th element in the continued correlation function of the k-th subarray, which consists of N-k+1 elements, and ${\mathit{w}}^{\left(k\right)}={\left[w_1^{\left(k\right)},w_2^{\left(k\right)},\dots ,w_{N-k+1}^{\left(k\right)}\right]}^T$denotes the noise eigenvector of the k-th subarray.

We now derive the relationship between the autoregressive coefficients $w_i^{\left(k\right)}$ in (86) for the Pisarenko and MUSIC methods under the assumption of a repeated noise eigenvalue.

The initial values of the correlation sequences in (86) are known and identical. Because the diagonal elements are not used in the prediction, these values can be drawn from either the covariance matrix $\mathbf{R}$ or its noise-free counterpart $\mathbf{C}=\mathbf{R}-{\sigma }^2\mathbf{I}$.

According to the singular extension theorem for Toeplitz covariance matrices [11], an N×N singular Toeplitz matrix C of rank d < N has a unique singular extension, provided its principal minor of order d is non-zero. The recurrence formula (86) specifies $\left(N-d\right)$ autoregressive extensions of the correlation function and their corresponding singular extensions. To verify the singularity of these extensions, we evaluate one step of the recurrence process for the autoregressive coefficient vector w_k. By assumption, w_k is the eigenvector corresponding to the zero eigenvalue of the matrix ${\mathbf{C}}^{\left(k\right)}={\mathbf{R}}^{\left(\mathit{k}\right)}-{\sigma }^2\mathbf{I}$:

$$\left[\begin{array}{cccc}{c}_1& c_2& \dots & c_{N-k+1}\\ c_2^{*}& c_1& \dots & c_{N-k}\\ \dots & \dots & \dots & \dots \\ c_{N-k+1}^{*}& c_{N-k}^{*}& \dots & c_1\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ \dots \\ w_{N-k+1}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \dots \\ 0\end{array}\right].$$

(87)

The extended matrix will be singular if its next extension—$\left(N-k+2\right)$-th column—is a linear combination of the preceding ones:

$$\left[\begin{array}{cccc}{c}_2& \dots & c_{N-k+1}& c_{N-k+2}\\ c_1& \dots & c_{N-k}& c_{N-k+1}\\ \dots & \dots & \dots & \dots \\ c_{N-k+1}^{*}& \dots & c_1& c_2\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ \dots \\ w_{N-k+1}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \dots \\ 0\end{array}\right].$$

(88)

The orthogonality of vector w_k to rows 2 through $N-k+1$ of the matrix in (88) follows directly from (87). Its orthogonality to the first row is a consequence of the recurrence formula (86) for the predicted matrix element. Consequently, each eigenvector (the vector of autoregression coefficients) w_k determines a singular extension of the noise-free covariance matrix for the $\left(N-k+1\right)$-element subarray.

While there exist $N-d$ infinite singular Toeplitz matrices of rank d that share a (d+1)×(d+1) block ${\mathbf{C}}^{(d+1)}$, they represent different extensions of the same submatrix ${\mathbf{C}}^{(d+1)}$. The singular extension theorem, however, guarantees uniqueness of the extension under these conditions. Therefore, all resulting autoregressive correlation sequences and their spectra are identical.

Let us now formulate the relationship between the Pisarenko and MUSIC spatial spectra. We define the spectrum of the k-th correlation sequence, which is generated by the noise eigenvector w_k in accordance with (93), as

$$P_{PSR,k}\left(\theta \right)={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{w}}_k\right|}^2}}.$$

(89)

Then (89) defining the MUSIC spatial spectrum estimate takes the following form:

$$P_{MUSIC}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^{N-d}}\frac{1}{P_{PSR,k}\left(\theta \right)}}}.$$

(90)

Thus, the MUSIC spectrum can be viewed as the harmonic mean of (N-d) Pisarenko spectra, while the Capon spectrum similarly equates to the harmonic mean of the maximum entropy estimates.

Since the subarray Pisarenko spectra of a uniform antenna array are identical, their average in the MUSIC estimate offers no advantage over a single Pisarenko spectrum. Thus, MUSIC and the Pisarenko method are identical when evaluating a single multiple noise eigenvalue. In the presence of a single multiple noise eigenvalue, the MUSIC spectral estimate becomes mathematically equivalent to the Pisarenko estimate.

While the MUSIC and Pisarenko spectra generally differ when noise eigenvalues vary or the antenna array is non-equidistant, the relationship in (90) holds true regardless of array geometry. However, for non-equidistant arrays, these spectra can no longer be interpreted as minimum entropy extensions of the correlation function.

The relationship between the MUSIC spectrum estimates and the minimum-entropy prediction of the correlation function is schematically illustrated in Figure 9.

Figure 9 compares the Pisarenko and MUSIC estimates from a conceptual standpoint, rather than focusing on specific computational scheme. To produce the Pisarenko spectral estimates, an equivalent inverse (non-recursive) filter with matching coefficients is used. Figure 10 illustrates two equivalent definitions of the Pisarenko spectrum.

Thus, the Capon and MUSIC methods represent harmonically averaged estimates of the spectra of infinite extensions of the correlation function beyond the observed region, corresponding to maximum entropy and minimum entropy (Pisarenko) extensions, respectively. For MUSIC and Pisarenko, these extensions equate to predicting the plane wavefronts of sources beyond the physical antenna aperture.

6. A Unified Framework for Super-Resolution Methods

6.1. Signal Models Spanning the Antenna Array and Beyond

In previous sections, we interpreted angular super-resolution methods as extrapolating received signal correlation functions or source wavefronts beyond the physical antenna aperture. To facilitate a comparison, Figure 11 provides a schematic summary of this concept. It illustrates signal models for the traditional and super-resolution methods, depicting linear antenna arrays and wavefronts both within and beyond the aperture for a single signal source.

Classical angular measurement methods utilize a plane-wave model for signal sources within the antenna aperture and zero elsewhere as presented in Figure 11(a). Variations of this approach incorporate spatial smoothing windows, which taper the signal amplitude transition at the aperture edges. This technique suppresses false peaks and sidelobe levels, reducing angular resolution.

Super-resolution algorithms that surpass the Rayleigh limit differ primarily in how they model wavefronts beyond the antenna aperture. Regardless of the specific technique, the ultimate performance of any super-resolution method depends on two key factors: how it extrapolates the correlation function outside the aperture, and how closely its underlying assumptions match the physical environment.

The maximum entropy method, illustrated in Figure 11(b), extrapolates the source wavefronts beyond the antenna aperture maximizing the entropy of their correlation function. Conversely, the Pisarenko method minimizes this entropy, thereby extending plane wavefronts beyond the aperture as shown in Figure 11(c).

Figure 11(d) shows the Capon signal model that averages N maximum entropy estimates of the source wavefronts, using antenna subarrays ranging from 1 to N elements.

Figure 11(e) shows the interpretation of the MUSIC method as extending plane signal wavefronts beyond the physical array. In the presence of d signal sources, N-d correlation functions are extrapolated using the minimum entropy criterion. Each extrapolated function corresponds to extending the plane wavefronts of all d signal sources. The resulting MUSIC spatial spectrum is the harmonic average of the spectra derived from the extended correlation functions.

Figure 11(f) illustrates a maximum entropy variant of ESPRIT, which we call Approximate ESPRIT. This variant retains the classical ESPRIT array configuration but utilizes a maximum entropy extension of the correlation function.

Operating on two identical subarrays, the ESPRIT method can be interpreted as extrapolating the plane wavefronts of signal sources beyond the physical array, as shown in Figure 11(g). This wavefront extrapolation is equivalent to applying a minimum entropy criterion to the extended blocks of the correlation function. The subarrays are equidistantly spaced. For conceptual clarity, they are depicted as non-overlapping, though they may overlap.

Figure 11(h) illustrates the linear prediction interpretation of the VESPA algorithm. In the first stage, fourth-order cumulants are used to synthesize a virtual secondary antenna array. The second stage then extrapolates the plain signal wavefronts to subsequent predicted virtual subarrays, mirroring the standard ESPRIT approach.

6.2. Common Space for Angular Super-Resolution Methods

We formulate angular super-resolution methods using a unified framework based on linear prediction and correlation-function extrapolation beyond the observed domain. This approach defines a multidimensional space of super-resolution estimates, whose axes represent entropy, averaging of estimates, and multichannel subarray prediction.

The origin of the estimation space (Figure 12) represents the Pisarenko method. This point features a minimum-entropy extrapolated correlation function, single-channel prediction, and no averaging of the estimates.

The extreme opposite of the entropy axis represents the maximum entropy method for single-channel, unaveraged prediction. Estimates with intermediate entropy values lie between the Pisarenko and maximum entropy limits.

The harmonic averaging of correlation function predictions corresponds to the MUSIC estimate for minimum-entropy predictions and the Capon estimate for maximum-entropy predictions. Specifically, the number of averaged predictions equals the number of signal sources ($N-d)$ for the MUSIC method, and the number of antenna array channels (N) for the Capon method. However, our graphical analysis indicates the presence or absence of averaging, without differentiating between the specific positions of these methods along the "Number of averaged predictions" axis.

The upper face of the cube represents the space of multichannel prediction by subarrays. The boundary methods in this plane are multidimensional counterparts to traditional single-channel techniques; for example, the ESPRIT algorithm corresponds to Pisarenko minimum entropy method, while approximate ESPRIT parallels the maximum entropy method. Points 1 and 2 designate the boundary states that have not yet been addressed by existing super-resolution algorithms. Details on deriving estimate 1 are provided below.

Our representation of the Pisarenko, maximum entropy, MUSIC, and Capon methods in the single-channel forecast plane is conceptually similar to the approach proposed in [27], which utilizes the signal-to-noise ratio rather than the entropy of the predicted correlation function. More importantly, however, the proposed angular super-resolution space not only establishes the relationships among these key method groups but also enables smooth transitions between them. This allows for the generation of arbitrary intermediate estimates that lie between established techniques.

Any point within the cubic region can be assigned an estimate, making this domain the estimate space for angular super-resolution methods. In Section VII we show how to practically obtain intermediate estimates bounded by this space. Furthermore, we introduce specific algorithms that enable for seamless transitions between estimates within a unified computational framework.

We now extend our super-resolution estimate space to higher-order statistics, requiring a redefinition of entropy. Since the standard log-determinant of a covariance matrix is valid only for Gaussian processes, evaluating non-Gaussian signals requires returning to the fundamental probability density definition. Without specifying this density, entropy cannot be directly expressed in terms of cumulant matrices. However, for extended signals, we are not interested in the exact entropy, but rather in how it reflects the nature of the spatial spectrum. Therefore, we introduce a surrogate measure analogous to the Gaussian covariance matrix determinant.

The cumulant matrix, used in the VESPA method to estimate the angular coordinates of sources, is the sum of the normalized source covariance matrices multiplied by specific scaling factors. In other words, the VESPA method leverages fourth-order statistics to calculate second-order statistics, which are then used to estimate the spectra or angular coordinates of the signal sources. In the presence of non-Gaussian uncorrelated noise, this cumulant matrix takes the form:

$$\tilde{\mathbf{R}}={\displaystyle \sum _{i=1}^d}{\beta }_i{\mathbf{R}}_i=\mathbf{A}\mathsf{\Gamma }{\mathbf{A}}^H+{\sigma }^2\mathbf{I},$$

(91)

where R_i is the covariance matrix of the i-th source signal across the antenna aperture; ${\beta }_i={\gamma }_i/{\sigma }_i^2$ is a scaling factor defined as the ratio of the fourth-order cumulant to the signal variance; ${\gamma }_i$ and ${\sigma }_i^2$ are the fourth-order cumulant and variance of the i-th source signal, respectively; Γ = diag{γ_i} is the diagonal matrix of these cumulants; and σ² is the variance of the uncorrelated non-Gaussian noise in the receiving channels.

We use the determinant of the cumulant matrix $\tilde{\mathbf{R}}$ as the third dimension in our space of angular super-resolution estimates based on fourth-order statistics.

As discussed in Section IV, the VESPA method performs multichannel prediction utilizing the minimum cumulant matrix determinant criterion without averaging of the estimates. Figure 13 illustrates the overall space of fourth-order super-resolution estimates, within which the VESPA method represents the only filled boundary cell. The remaining vertices bounding the cumulant estimate space are defined by the natural physical limits of the key parameters: the full cumulant matrix determinant, the number of simultaneously predicted channels, and the number of averaged predictions.

In the single-channel prediction plane, Figure 13 illustrates boundary estimates for non-Gaussian signals, mirroring the second-order statistics approach shown in Figure 12. Clearly, using a cumulant matrix instead of a covariance matrix allows for the development of higher-order counterparts to standard correlation-based methods.

6.3. An Example of Synthesizing the ‘Missing’ Super-Resolution Method

Some cube vertices in the angular super-resolution space (Figure 12) are empty. To show that a corresponding super-resolution method can be derived for every vertex or internal point, we get an estimate for the method labeled 1.

Spatial spectrum estimates for unfilled cube vertices can be synthesized via adjacent edge methods. For vertex 1, you can use MUSIC (with a multichannel version) or ESPRIT (with harmonic averaging). We focus on the ESPRIT approach.

Standard ESPRIT uses only two identical subarrays. Applying harmonic averaging, however, requires more than two subarrays. We therefore use K identical physical subarrays, labeled as AР₁…AР_К in Figure 14 (which illustrates the signal models on and beyond the antenna aperture).

Without source cross-correlation, the real antenna array's covariance matrix for ESPRIT becomes:

$$\mathbf{C}=\left[\begin{array}{cccc}{\mathbf{C}}_1& {\mathbf{C}}_2& \dots & {\mathbf{C}}_k\\ {\mathbf{C}}_1^H& {\mathbf{C}}_1& \dots & {\mathbf{C}}_{k-1}\\ \dots & \dots & \dots & \dots \\ {\mathbf{C}}_k^H& {\mathbf{C}}_{k-1}^H& \dots & {\mathbf{C}}_1\end{array}\right],$$

(92)

where C_i is the mutual covariance matrix of signals of the m-th and $\left(m+i-1\right)$-th subarrays; K is the number of real subarrays.

We assume the number of signal sources d satisfies $d\le N$, where N is the dimension of matrix C_i (the number of channels in a subarray). Following Section III, the generalized eigenvalues of matrix pencils formed by adjacent matrix pairs C_i and C_i+1 matrix C are identical. These eigenvalues directly yield the angular directions of the signal sources:

$${\mathbf{C}}_1\mathbf{U}={\mathbf{C}}_{i+1}\mathbf{U}\mathsf{\Lambda }$$

(93)

where U = [u₁, u₂, ..., u_N], and $\mathsf{\Lambda }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}{\left\{{\lambda }_m\right\}}_{m=1,\dots ,N}$, u_m, λ_m represents the m-th generalized eigenvector and the m-th generalized eigenvalue of the beam, respectively. The angular coordinates of the sources relate to the generalized eigenvalues as follows:

$$\mathsf{\Lambda }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\exp \left(j2\pi \left|\mathsf{\Delta }\right|\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_1\right),\dots ,\exp \left(j2\pi \left|\mathsf{\Delta }\right|\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_d\right),0,\dots ,0\right\}.$$

(94)

When matrices U and Λ satisfy equation (93), they can be used to construct a matrix of eigenvectors for C corresponding to the zero eigenvalue. The block columns of this matrix, each containing K blocks of size $N\times N$, can then be interpreted as the block eigenvectors of the noise-free block covariance matrix C:

$$\mathrm{C}\mathrm{B}=0$$

(95)

or

$$\left[\begin{array}{cccccc}{C}_1& C_2& C_3& \dots & C_{K-1}& C_K\\ C_2^H& C_1& C_2& \dots & C_{K-2}& C_{K-1}\\ C_3^H& C_2^H& C_1& \dots & C_{K-3}& C_{K-2}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ C_{K-1}^H& C_{K-2}^H& C_{\mathrm{K}-3}^H& \dots & C_1& C_2\\ C_K^H& C_{K-1}^H& C_{K-2}^H& \dots & C_2^H& C_1\end{array}\right] \left[\begin{array}{cccc}{B}_1& 0& \dots & 0\\ B_2& B_1& \dots & 0\\ 0& B_2& \dots & 0\\ \dots & \dots & \dots & \dots \\ 0& 0& \dots & B_1\\ 0& 0& \dots & B_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ \dots \\ 0\\ 0\end{array}\right],$$

where B₁ and B₂ are $N\times N$ square matrices. Since these matrices are not uniquely determined, one valid choice satisfying equation (93) is:

 

(96)

Each block eigenvector forms a multichannel prediction vector. When the full antenna array has a block-Toeplitz covariance matrix, these predictions match. If cross-correlated sources break this block-Toeplitz structure, the predictions diverge, changing the structure of matrix B.

Therefore, the spectrum estimate for the proposed method is

$$P_1\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^{K-1}}{‖{\mathbf{a}\left(\theta \right)}^H{\mathbf{B}}^{\left(k\right)}‖}^2}},$$

(97)

where B^(k) is the k-th block column of the matrix B according to (95).

Clearly, if matrices B₁ and B₂ are defined by (96), the spectral peaks in (97) correspond to the generalized eigenvalues of the matrix pencil $\left\{{\mathbf{C}}_i,{\mathbf{C}}_{i+1}\right\}$ because

$$P_1\left(\theta \right)={\displaystyle \frac{1}{{‖{{\mathbf{a}}_1\left(\theta \right)}^H\left[{\mathbf{B}}_1-z\left(\theta \right){\mathbf{B}}_2\right]‖}^2}}={\displaystyle \frac{1}{{‖{{\mathbf{a}}_1\left(\theta \right)}^H\left[\mathbf{U}-z\left(\theta \right)\mathbf{U}\mathsf{\Lambda }\right]‖}^2}},$$

(98)

where $z\left(\theta \right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-j2\pi \left|∆\right|\mathrm{s}\mathrm{i}\mathrm{n}\theta \right)$; a₁(θ) is the steering vector at the aperture of the first subarray.

On the other hand, the estimate in (97) represents the harmonic average $\left(K-1\right)$ of the ESPRIT estimates:

$$P_1\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^{K-1}}\frac{1}{P_{ESP,k}\left(\theta \right)}}},$$

(99)

where

$$P_{ESP,k}\left(\theta \right)={\displaystyle \frac{1}{{‖{\mathbf{a}\left(\theta \right)}^H{\mathbf{B}}^{\left(k\right)}‖}^2}}.$$

(100)

When the array's covariance matrix is block-Toeplitz, estimates (100) from all eigenvector blocks B^(k) coincide. In contrast, non-Toeplitz blocks are unique, meaning harmonic averaging (99) produces an estimate different from standard ESPRIT.

Thus, estimate (97) is a multichannel prediction estimate with minimal entropy (extension of plane wavefronts beyond the antenna aperture) and harmonic averaging of predictions. This makes it a multichannel version of the MUSIC method.

7. Unifying Super-Resolution Implementations

Placing super-resolution methods in a common space lets us parameterize spatial spectra. As a result, we can obtain any single-channel prediction (the bottom plane in Figure 12) through one universal linear filtering process.

From a computational view, we can estimate spatial spectrum parameters in two ways: using the antenna array signal covariance matrix or deriving the parameters directly from the channel signal vectors. Every orthogonalization procedure for sample snapshot vectors directly corresponds to a specific transformation of the covariance matrix. For instance, the Gram-Schmidt method maps to Gaussian elimination, and the modified Gram-Schmidt method maps to Cholesky factorization of a Hermitian matrix. As a baseline, we consider orthogonalization procedures applied to the signal snapshots, while also demonstrating equivalent procedures based on covariance matrices.

7.1. J-Orthogonalization

Spatial spectra can be computed using either covariance matrix operations or direct processing of received signals. The second approach derives spectra through sequential orthogonalization of the antenna array signals, such as Gram-Schmidt orthogonalization, lattice filters, or unitary rotations [15,29,30,31]. For clarity, we will demonstrate the J-orthogonalization principle using a modified Gram-Schmidt algorithm [32].

Since the objects to be orthogonalized are vectors of the signal snapshots in the antenna array channels, we redefine the sample covariance matrix for a block of samples. We use the maximum likelihood estimate, omitting the scaling factor as it has no impact on angular coordinate estimation:

$$\mathbf{R}=\mathbf{X}{\mathbf{X}}^H,$$

(101)

where X = [X₁, X₂, ..., X_N]^T is the matrix of signal samples across the antenna array channels; X_n is the K-dimensional vector representing the complex signal amplitudes at the output of the n-th channel, where K is the sample size; and N is the number of channels in the array.

Let Q be a non-singular square matrix that linearly transforms the original signals X₁, X₂, ..., X_N into an orthonormal basis Z₁, Z₂, ..., Z_N,

$$\mathbf{Z}=\mathbf{Q}\mathbf{X};\mathbf{Z}{\mathbf{Z}}^H=\mathbf{I}.$$

(102)

In this case, the matrix Q is the square root of the inverse covariance matrix R^-1,

$$\mathbf{Q}\mathbf{X}{\mathbf{X}}^{\mathit{H}}{\mathbf{Q}}^H=\mathbf{I};\mathbf{X}{\mathbf{X}}^{\mathit{H}}={\left(\mathbf{Q}{\mathbf{Q}}^H\right)}^{-1};{\mathbf{R}}^{-1}=\mathbf{Q}{\mathbf{Q}}^H.$$

(103)

Without restrictions on the Q matrix, multiple decorrelation transformations satisfy this condition, including the modified Gram-Schmidt orthogonalization algorithm. In this case, the impulse response of the tuned orthogonalizing filter is the lower triangular square root of the inverse covariance matrix,

$${\mathbf{R}}^{-1}=\mathbf{L}{\mathbf{L}}^H,$$

(104)

representing a specific version of the matrix Q.

The Gram-Schmidt orthogonalization filter yields super-resolution estimates represented by the columns or rows of the lower triangular matrix L, including the maximum entropy and Capon estimates.

Next, we introduce the J-orthogonalization procedure, a generalization of the standard orthogonalization. We will use this approach to derive maximum and minimum entropy spatial spectrum estimates.

We call two vectors x and y J-orthogonal if

$${\mathbf{x}}^H\mathbf{J}\mathbf{y}=\mathbf{0},$$

(105)

where J is a square matrix.

We supplement the received signal snapshot matrix X with a scaled identity matrix $\widehat{\sigma }\mathbf{I}$, defining the new matrix as M,

$$\mathbf{M}=\left[\begin{array}{cc}\mathbf{X}& \widehat{\sigma }\mathbf{I}\end{array}\right].$$

(106)

Instead of orthogonalizing the rows of X as in (102), we enforce J-orthogonality on the rows of the matrix M,

$$\mathbf{Z}=\mathbf{Q}\mathbf{M};\mathbf{Z}\mathbf{J}{\mathbf{Z}}^H=\mathbf{I},$$

(107)

where Q denotes the transformation matrix, and

 

(108)

where I_K and I_N denote the identity matrices of sizes K×K and N×N, respectively.

This determines the set of vectors (106), the J-orthogonalization condition (107), and the basis J (108). Like standard orthogonalization, various J-orthogonalization procedures exist. While mathematically equivalent, they differ in numerical stability and some other features. We apply the modified Gram-Schmidt algorithm to orthogonalize the rows of M in the basis J as follows:

$${\mathbf{M}}^{(n+1)}=\left[\begin{array}{cc}{\mathbf{X}}^{(n+1)}& {\mathbf{H}}^{(n+1)}\end{array}\right]={\tilde{\mathbf{F}}}_n{\mathbf{M}}^{\left(n\right)}={\tilde{\mathbf{F}}}_n...{\tilde{\mathbf{F}}}_1\mathbf{M};n=1,\dots ,N-1,$$

(109)

$${\tilde{\mathbf{F}}}_n=\left[\left(\begin{array}{ccccccc}1& & & & & & 0\\ & \ddots & & & & & \\ & & 1& & & & \\ & & & \left(1/\sqrt{{\tilde{g}}_n}\right)& & & \\ & & & {\tilde{f}}_{n+1,n}& 1& & \\ & & & ⋮& & \ddots & \\ 0& & & {\tilde{f}}_{N,n}& & & 1\end{array}\right)\right],$$

$${\tilde{g}}_n={\displaystyle \sum _{k=1}^K}{\left|x_{n,k}^{\left(n\right)}\right|}^2-{\displaystyle \sum _{k=1}^K}{\left|h_{n,k}^{\left(n\right)}\right|}^2={‖{\mathbf{X}}_n^{\left(n\right)}‖}^2-{‖{\mathbf{H}}_n^{\left(n\right)}‖}^2,$$

(110)

$${\tilde{f}}_{i,n}=-{\displaystyle \frac{{\displaystyle \sum _{k=1}^K}{x}_{n,k}^{\left(n\right)*}{x}_{i,k}^{\left(n\right)}-{\displaystyle \sum _{k=1}^K}{h}_{n,k}^{\left(n\right)*}{h}_{i,k}^{\left(n\right)}}{{\displaystyle \sum _{k=1}^K}{\left|x_{n,k}^{\left(n\right)}\right|}^2-{\displaystyle \sum _{k=1}^K}{\left|h_{n,k}^{\left(n\right)}\right|}^2}}={\displaystyle \frac{{\mathbf{X}}_n^{\left(n\right)H}{\mathbf{X}}_i^{\left(n\right)}-{\mathbf{H}}_n^{\left(n\right)H}{\mathbf{H}}_i^{\left(n\right)}}{{‖{\mathbf{X}}_n^{\left(n\right)}‖}^2-{‖{\mathbf{H}}_n^{\left(n\right)}‖}^2}},i=n+1,\dots ,N.$$

(111)

Applying the J-orthogonalization procedure yields J-orthogonal rows for M at the filter outputs (Figure 15):

M^(N)JM^(N)H = I.

(112)

The J-orthogonalization filter's impulse response is still represented by a lower triangular matrix:

$$\tilde{\mathbf{L}}=\tilde{\mathbf{D}}{\tilde{\mathbf{F}}}_{N-1}\dots {\tilde{\mathbf{F}}}_2{\tilde{\mathbf{F}}}_1,$$

(113)

where $\tilde{\mathbf{D}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{\displaystyle \frac{1}{{\tilde{g}}_1}},\dots ,{\displaystyle \frac{1}{{\tilde{g}}_N}}\right\}$.

Suppose, the antiregularization parameter ${\widehat{\sigma }}^2$ differs from all eigenvalues of the sample covariance matrix, ensuring the matrix is non-singular. Then, condition (107) for J-orthogonality can be expressed as:

$$\tilde{\mathbf{L}}\mathbf{M}\mathbf{J}{\mathbf{M}}^H{\tilde{\mathbf{L}}}^H=\tilde{\mathbf{L}}(\mathbf{X}{\mathbf{X}}^H-{\sigma }^2\mathbf{I}{)\tilde{\mathbf{L}}}^H=\mathbf{I}.$$

(114)

Thus, the matrix is not the square root of the inverse covariance matrix R^-1, but rather of the matrix ${\mathbf{R}}_{\sigma }^{-1}$,

$${\mathbf{R}}_{\sigma }^{-1}={\left(\mathbf{X}{\mathbf{X}}^H-{\widehat{\sigma }}^2\mathbf{I}\right)}^{-1}={\tilde{\mathbf{L}}}^H\tilde{\mathbf{L}}.$$

(115)

Here, ${\widehat{\sigma }}^2$acts as an anti-regularization parameter. While regularization adds a value ${\widehat{\sigma }}^2$ to the diagonal (${\mathbf{R}}_{\sigma }=\mathbf{X}{\mathbf{X}}^H+{\widehat{\sigma }}^2\mathbf{I}$) to improve matrix conditioning, anti-regularization subtracts it (${\mathbf{R}}_{\sigma }=\mathbf{X}{\mathbf{X}}^H-{\widehat{\sigma }}^2\mathbf{I}$) to do the opposite. The extreme case where σ² equals the noise eigenvalue is also important.

The Gram-Schmidt J-orthogonalization process (111) can also be applied in a regularization mode by changing the subtraction in both the numerator and denominator to addition.

Knowing when J-orthogonalization is feasible is essential. While standard orthogonalization needs the sample covariance matrix R being full rank, J-orthogonalization strictly requires that neither R nor its leading principal submatrices eigenvalues equal to ${\widehat{\sigma }}^2$. However, setting ${\widehat{\sigma }}^2$ to precisely match a noise eigenvalue of R allows J-orthogonalization to generate the eigenvectors needed in projection (minimum-entropy) methods such as Pisarenko and MUSIC.

When the parameter ${\widehat{\sigma }}^2$ is set to the minimum (noise) eigenvalue of the sample covariance matrix ${\widehat{\sigma }}^2={\lambda }_{min}$, then the rank of ${\mathbf{R}}_{\sigma }$ corresponds to the number of signal sources, d,

 

(116)

We begin with two lemmas on the normalizing factors under condition (116) before proving the theorem.

Lemma 1.

For d = rank{R_σ}, the factor g_d+1 vanishes at step (d+1) of the J-orthogonalization process (109)–(111),

$${g_{d+1}=‖{\mathbf{X}}_{d+1}^{(d+1)}‖}^2-{‖{\mathbf{H}}_{d+1}^{\left(d+1\right)}‖}^2=0.$$

(117)

Proof. 

Assume for contradiction that ${‖{\mathbf{X}}_{d+1}^{(d+1)}‖}^2-{‖{\mathbf{H}}_{d+1}^{(d+1)}‖}^2\ne 0$. This implies that the leading submatrix of order (d+1), given by ${\mathbf{R}}_{\sigma }^{(d+1)}={\mathbf{R}}^{(d+1)}-{\widehat{\sigma }}^2{\mathbf{I}}_{d+1}$, is of full rank. Consequently, rank{R_σ} ≥ d+1. This directly contradicts the initial assumption that R_σ = d, where d is the number of signal sources.

Lemma 2. 

The normalizing factor g_n is nonzero for the first d steps of the J-orthogonalization process,

$${g_n=‖{\mathbf{X}}_n^{\left(n\right)}‖}^2-{‖{\mathbf{H}}_n^{\left(n\right)}‖}^2\ne 0,n=1,\dots ,d.$$

(118)

Proof. 

Assume that ${‖{\mathbf{X}}_n^{\left(n\right)}‖}^2-{‖{\mathbf{H}}_n^{\left(n\right)}‖}^2=0$ for all n = 1, ..., d. This implies that the noise-free covariance matrix of the subarray containing the first n antenna elements, ${\mathbf{R}}_{\sigma }^{\left(n\right)}={\mathbf{R}}^{\left(n\right)}-{\widehat{\sigma }}^2{\mathbf{I}}_n$, is singular. However, for n < d, this matrix is of full rank, regardless of the noise component, which leads to a contradiction.

Theorem. 

Let λ_min be the minimum eigenvalue of the matrix R = XX^H. The vectors ${\mathbf{H}}_j^{\left(d+1\right)}$ (j = d+1, ..., N) obtained by the J-orthogonalization (109)–(111) of the rows of the matrix $\mathbf{M}=\left[\begin{array}{cc}\mathbf{X}& \widehat{\sigma }\mathbf{I}\end{array}\right]$ with parameter $\widehat{\sigma }=\sqrt{{\lambda }_{\mathrm{m}\mathrm{i}\mathrm{n}}}$ are the eigenvectors of R corresponding to λ_min.

Proof.

The J-orthogonalization process (109)–(111) acts as a linear filtering procedure for the vectors in packet M. At the filter input, the square block ${\mathbf{H}}^1$ of matrix M is a scalar matrix ${\mathbf{H}}^1=\sigma \mathbf{I}$. After the n-th step, it becomes a transformation matrix representing the filter's impulse response—up to the coefficient $\widehat{\sigma }$. Consequently, vectors ${\mathbf{X}}_j^{\left(d+1\right)}$ equal the product of the input signal packet X and the corresponding impulse response,

$${\mathbf{X}}_j^{\left(d+1\right)H}{\mathbf{X}}_j^{\left(d+1\right)}={\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{H}}_j^{\left(d+1\right)},$$

(119)

where ${\mathit{e}}_k$ is the k-th noise eigenvector of the covariance matrix and d is the number of signal sources.

Passing steering vectors through a tuned Gram-Schmidt orthogonalization filter yields maximum entropy or Capon spatial spectrum estimates. These estimates come from inverting the output power of the N-th filter output (124) for maximum entropy, or the sum of all filter outputs (125) for Capon.

Another method directly obtains the desired l_*k vectors in the J-orthogonalization filter with parameter $\widehat{\sigma }=1$, adjusting the filter coefficients according to the standard Gram-Schmidt orthogonalization algorithm. The 'tails' of the packet M are excluded from the coefficient calculations; instead, they are filtered along with the packet's training portion and then substituted into the maximum entropy (124) and Capon (125) spectrum estimates correspondingly:

$${\mathbf{X}}_j^{\left(d+1\right)H}{\mathbf{X}}_j^{\left(d+1\right)}={\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{H}}_j^{\left(d+1\right)}.$$

(120)

Using (119) in (120), we find

$${\displaystyle \frac{1}{\sigma }}{\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{X}{\mathbf{X}}^H\mathbf{H}}_j^{\left(d+1\right)}={\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{H}}_j^{\left(d+1\right)},{\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{R}\mathbf{H}}_j^{\left(d+1\right)}={\sigma }^2{‖{\mathbf{H}}_j^{\left(d+1\right)}‖}^2,$$

(121)

$${\displaystyle \frac{{\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{R}\mathbf{H}}_j^{\left(d+1\right)}}{{‖{\mathbf{H}}_j^{\left(d+1\right)}‖}^2}}={\widehat{\sigma }}^2.$$

(122)

The left side of the last equation is the Rayleigh ratio. Its minimum value equals the smallest eigenvalue of the numerator's Hermitian matrix,

$$\underset{{\mathbf{H}}_j^{\left(d+1\right)}}{\min }\left\{{\displaystyle \frac{{\mathbf{H}}_j^{\left(d+1\right)H}{\mathbf{R}\mathbf{H}}_j^{\left(d+1\right)}}{{‖{\mathbf{H}}_j^{\left(d+1\right)}‖}^2}}\right\}={\lambda }_{min}.$$

(123)

As ${\widehat{\sigma }}^2$ in (122) is the minimum eigenvalue of R by assumption, it follows that ${\mathbf{H}}_j^{\left(d+1\right)}$ is the associated eigenvector, completing the proof.

Thus, the vectors ${\mathbf{H}}_j^{\left(d+1\right)}$ (j = d+1, ..., N) represent the noise eigenvectors of the covariance matrix, assuming the parameter $\widehat{\sigma }$ is the square root of the minimum noise eigenvalue. This enables the implementation of the Pisarenko and MUSIC minimum entropy methods within the J-orthogonalization process.

7.2. Bridging Minimum and Maximum Entropy Spatial Spectra

We show that a J-orthogonalization filter yields the angular super-resolution estimates for all single-channel prediction methods on the lower bound of the estimate space cube (Figure 12). These methods include Pisarenko, MUSIC, Capon, and maximum entropy. First, we rewrite each spatial spectrum expression to fit an adaptive linear filter. We omit normalizing factors since they do not change the DOA estimates, unless doing so causes confusion.

The maximum entropy spatial spectrum estimate is defined as:

$$P_{ME}\left(\theta \right)={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathit{l}}_{*1}\right|}^2}},$$

(124)

where a(θ) is the steering vector, and l_*1 is the first column of the lower triangular matrix L from the decomposition R^-1 = LL^H.

The Capon estimate is given by:

$$P_{CP}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle \sum _{k=1}^N}{\left|{\mathbf{a}\left(\theta \right)}^H{\mathit{l}}_{*k}\right|}^2}},$$

(125)

where l_*k is the k-th column of the matrix L

The Pisarenko spatial spectrum estimate is:

$$P_{PSR}\left(\theta \right)={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{e}}_1\right|}^2}},$$

(126)

where e₁ is the unique noise eigenvector of the covariance matrix.

The MUSIC estimate is:

$$P_{MUSIC}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle \sum _{k=1}^{N-d}}{\left|{\mathbf{a}\left(\theta \right)}^H{\mathit{e}}_k\right|}^2}}.$$

(127)

where ${\mathit{e}}_k$ is the k-th noise eigenvector of the covariance matrix and d is the number of signal sources.

Passing steering vectors through a tuned Gram-Schmidt orthogonalization filter yields maximum entropy or Capon spatial spectrum estimates. These estimates come from inverting the output power of the N-th filter output (124) for maximum entropy, or the sum of all filter outputs (125) for Capon.

Another method directly obtains the desired l_*k vectors in the J-orthogonalization filter with parameter $\widehat{\sigma }=1$, adjusting the filter coefficients according to the standard Gram-Schmidt orthogonalization algorithm. The 'tails' of the packet M are excluded from the coefficient calculations; instead, they are filtered along with the packet's training portion and then substituted into the maximum entropy (124) and Capon (125) spectrum estimates correspondingly:

$$P_{ME}\left(\theta \right)={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{H}}_N^{\left(N\right)}\right|}^2}},$$

(128)

$$P_{CP}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^N}{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{H}}_k^{\left(N\right)}\right|}^2}}.$$

(129)

Next, we consider obtaining minimum entropy estimates based on the eigenvectors of the noise subspace during J-orthogonalization. To transition from maximum entropy and Capon estimates to the minimum entropy estimates of the Pisarenko and MUSIC, it is sufficient to set the antiregularization parameter ${\widehat{\sigma }}^2$ in the J-orthogonalization filter (109)–(111) equal to the minimum (noise) eigenvalue λ_min of the covariance matrix.

The precise equality of the antiregularization parameter to the noise eigenvalue is less critical than it seems. Choosing an underestimated value only shifts the spatial spectrum toward higher entropy. As shown below, the spectrum changes smoothly across the antiregularization parameter from λ_m to zero, moving from the minimum to the maximum entropy spectrum. The following estimate of the spatial spectrum using the outputs of the J-orthogonalization filter corresponds to the Pisarenko method:

$$P_{PSR}\left(\theta \right)={\displaystyle \frac{1}{{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{H}}_N^{(d+1)}\right|}^2}}.$$

(130)

The MUSIC method uses $N-d$ noise subspace vectors from the covariance matrix. We assume the noise eigenvalues are equal, meaning the minimum eigenvalue has a multiplicity of $N-d$. These noise vectors appear as $N-d$ vectors ${\mathbf{H}}_j^{(d+1)}$ (j = d+1, ..., N) at the output of the d-th filter cascade (Figure 15) using J-orthogonalization. Because they correspond to a multiple eigenvalue, these vectors are not mutually orthogonal yet. However, they fully define the noise subspace. We orthogonalize them across the remaining $N-d-1$ filter cascades using the standard Gram-Schmidt algorithm without changing the filter structure. This yields mutually orthogonal vectors defining the MUSIC spatial spectrum:

$$P_{MUSIC}\left(\theta \right)={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^{N-d}}{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{H}}_{(d+j)\perp }^{\left(N\right)}\right|}^2}}.$$

(131)

By applying the J-orthogonalization technique, we derive four primary super-resolution estimates situated at the opposite extremes of the entropy scale. We will next develop intermediate-entropy estimates that bridge the gap between the Pisarenko and maximum entropy spectrums, as well as and the MUSIC and Capon spectrums.

Producing spectrum estimates based on intermediate entropy values in the J-orthogonalization filter is equivalent to obtaining Capon estimates. However, the antiregularization parameter ${\widehat{\sigma }}^2$ must range between zero and the noise eigenvalue of the covariance matrix matrix $0<{\widehat{\sigma }}^2<{\lambda }_{\mathrm{m}\mathrm{i}\mathrm{n}}$.

Varying the number of averaged predictions gives us every super-resolution estimate in the single-channel prediction plane (Figure 12):

$$P\left(\theta \right){|}_{\sigma ,K}={\displaystyle \frac{1}{{\displaystyle {\sum }_{k=1}^K}{\left|{\mathbf{a}\left(\theta \right)}^H{\mathbf{H}}_{N-K+k}^{\left(N\right)}\right|}^2}},$$

(132)

where K denotes the number of averaged predictions.

We have demonstrated how to obtain all super-resolution estimates, from minimum to maximum entropy, including the optional harmonic averaging of predictions within the J-orthogonalization filter. Figure 16 illustrates the derivation of spectrum estimates for the four super-resolution methods on the single-channel prediction plane (Figure 12), including any intermediate estimates. In all scenarios, the filter inputs consist of a composite packet M, which contains signal samples from the antenna array channels X, and a scalar square matrix σI.

When the parameter σ equals one and the filter executes a standard Gram-Schmidt orthogonalization on the receiving signals X, filtering the entire augmented matrix M = [X | I], we obtain maximum entropy and Capon estimates.

Pisarenko and MUSIC estimates are derived by setting the filter's antiregularization parameter equal to the covariance matrix's minimum eigenvalue. As filtering proceeds, the last $N-d$ rows of the submatrix $\widehat{\sigma }\mathbf{I}$ become eigenvectors of the noise subspace.

Spatial spectrum estimates can be computed directly from the prediction parameter vectors of the J-orthogonalization filter. Alternatively, they can be obtained by passing the steering vectors a(θ_i) through the tuned filter, an approach illustrated in Figure 17.

We demonstrate how a J-orthogonalization filter can yield maximum entropy, Pisarenko, Capon, and MUSIC super-resolution estimates, and their intermediate variants. Our simulation features a 16-element ULA receiving signals from four uncorrelated Gaussian sources in additive uncorrelated Gaussian noise. The angular coordinates of these sources are set to -0.2, -0.16, 0.1, and 0.25, with corresponding powers of 10, 5, 10, and 1 dB.

Figure 18 illustrates the spatial spectra generated by the J-orthogonalization filter for the maximum entropy, Capon, Pisarenko, and MUSIC estimators.

Figure 19 and Figure 20 display the spatial spectrum estimates across a range of intermediate entropy values from maximum to minimum, without and with harmonic averaging, respectively. The entropy axis shows how the spatial spectrum changes across different values of the antiregularization parameter (σ² = 0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0). Position 7 in the foreground represents maximum entropy (σ² = 0), while position 1 in the background represents minimum entropy (σ² = 1). While the true covariance matrix has a minimum multiple noise eigenvalue of exactly 1, the sample covariance matrix yields varying noise eigenvalues ranging from 0.8 to 1.13. This means the Pisarenko and MUSIC methods use approximate extensions of the minimum-entropy spatial correlation function rather than exact extensions—which are equivalent to noise subspace projectors—their spatial spectra may differ slightly from Figure 19.

Figure 19 shows that the angular resolution increases as the antiregularization parameter approaches the noise eigenvalue (i.e., as entropy decreases toward Pisarenko estimate). This improvement is evident from the distinct resolution of the closely spaced sources at coordinates –0.2 and –0.16.

Figure 20 shows the spatial spectra obtained from the J-orthogonalization filter using harmonic averaging of the estimates. The Capon spectrum appears in the foreground, characterized by maximum entropy (zero antiregularization parameter, σ² = 0). Moving along the entropy axis, the final plot converges toward the MUSIC minimum entropy estimate (the antiregularization parameter σ² = 1).

8. Conclusions

Super-resolution methods can be organized within a coordinate space defined by three core factors: entropy, estimate averaging, and multichannel subarray prediction. This framework not only structures existing methods but also accommodates any intermediate estimates. Building on this, our approach leverages a universal J-orthogonalization algorithm to simultaneously generate spectral estimates across multiple super-resolution methods with virtually no extra computational overhead. Consequently, it enables the adaptive selection or combination of super-resolution methods to optimize performance under varying signal-to-noise ratios and array calibration conditions.