Low-Complexity Compressed Sensing-Aided Coherent Direction-of-Arrival Estimation for Large-Scale Lens Antenna Array

This paper delves into a novel compressed sensing (CS) strategy for estimating the directions of incoming signals in a coherent environment using a lens antenna array (LAA). In comparison to the well-known subspace-based algorithm family, CS techniques, such as the conventional orthogonal matching pursuit (COMP), can effectively address the direction-of-arrival (DoA) estimation problem requiring prior knowledge about the number of signals and offer lower complexity. However, they are susceptible to noise and can be adversely affected by multipath distortion. Leveraging the energy-concentrating property of an LAA, we first introduce the signal covariance matrix-based OMP (SCM-OMP) method that enhances the angular estimation performance, even in low-SNR regions. Subsequently, we propose the multiple sub-covariance matrices-based OMP (MSCM-OMP) to achieve a reduction in computational complexity. Simulation results demonstrate that the MSCM-OMP scheme also outper-forms other high-resolution DoA estimation methods.


I. INTRODUCTION
The lens antenna array (LAA) serves as an effective implementation of massive multiple-input-multiple-output (MIMO) technology, a pivotal enabler for millimeter-wave (mmWave) and terahertz (THz) wireless communication systems.Specifically, LAAs have garnered considerable attention due to their unique ability to manipulate electromagnetic (EM) waves, concentrate signal energy, achieve a high level of integration, and deliver superior energy efficiency compared to traditional discrete antenna arrays [1].Therefore, LAAs have the potential to be widely applied in mmWave and THz communications [1], precise tracking and localization tasks [2], and efficient wireless power transfer [3].It is imperative to underscore the fundamental role played by accurate direction-of-arrival (DoA) in these contexts.Precise DoA estimation empowers the selection of the discrete Fourier transform (DFT) beam, which, in turn, optimally aligns with the dominant signal path.This alignment proves indispensable for tasks such as spatial multiplexing, beamforming, interference mitigation, and energy-efficient communication.This significantly contributes to elevating signal quality, enhancing system capacity, and elevating the overall efficiency of the communication systems.
In recent years, various methodologies have emerged to address DoA estimation using an LAA [4]- [8].The subspacebased algorithm family, such as the conventional multiple signal classification (MUSIC) or estimation of signal parameters via the rotational invariance technique (ESPRIT), can achieve a high-resolution direction estimation performance [7].However, a common limitation of these approaches is their reliance on a priori knowledge of the number of incoming signals, necessitating additional computational resources for precise signal enumeration.It is important to underscore that an incorrect estimation of signal number results in significant errors in subspace-fitting-based DoA estimation.Furthermore, the conventional MUSIC encounters challenges when dealing with scenarios involving multipath signals.The spatial smoothing technique, namely the full-row Toeplitz matrices reconstruction (FTMR) [4], [9], [10], serves as an effective preprocessing tool to mitigate the issues posed by coherent signals.In addition, the authors in [8] propose an atomic norm minimization (ANM) approach to recover the original signal, then employ the MUSIC scheme to construct the spatial spectrum.It is noteworthy that these procedures, including ANM methods and other smoothing techniques, however, demand substantial computational resources.
Unlike traditional subspace-based methods, compressed sensing (CS) approaches tackle the problem of uncorrelated or coherent DoA estimation without modifying the data covariance matrix.Among the CS sparse signal recovery algorithms, the orthogonal matching pursuit (OMP) [11] stands out for its significantly lower computational complexity and faster processing speed [12].OMP-based DoA estimation focuses on reconstructing a spatial spectrum by leveraging the sparsity of non-zero values in the angular domain.However, despite its effectiveness in rapidly recovering DoA information, it is noted that the OMP-aided estimation process is not immune to some practical challenges.Specifically, the recovered signals may contain elements of noise, which can influence the precision of DoA estimates, especially in a low-SNR regime.In this paper, two frameworks inherited from the traditional CS technique are proposed to estimate the directions of incoming signals.First, the signal covariance matrix-based OMP algorithm, which utilizes the columns corresponding to active antennas, is proposed to improve the estimation performance in a highly contaminated environment.Next, the multiple sub-covariance matrices-based OMP (MSCM-OMP) is further introduced to reconstruct the partial spectra corresponding to different submatrices.Moreover, we mathematically analyze the MSCM-OMP scheme to prove its reduction in computational complexity as the size of the LAA or the number of antennas increases.The number of floating-point operations (FLOPs) [13] of MSCM-OMP is also derived.Last but not least, the simulation results demonstrate that our proposed MSCM-OMP not only improves estimation performance but also exhibits lower complexity compared to other state-of-theart algorithms.
Notations: The notations (•) T , (•) * , and (•) H indicate the transpose, conjugate, and conjugate transpose of a complex vector or matrix, respectively.The mth column or row of matrix A are denoted by A :,m or A m,: , respectively.Furthermore, a m represents the mth element of a vector a.The notation E {•} signifies the expectation of a random variable.In addition, the operator diag {•} with a square matrix inside generates a vector containing the diagonal elements of the input matrix.I M denotes an M ×M matrix with ones along the main diagonal and zeros elsewhere.Moreover, |A| indicates the cardinality of set A.

A. Coherent DoA Estimation Problem
An LAA is composed of a planar EM lens placed on the y-axis with physical length D y , and antenna elements placed on the focal arc.We denote m as the index of the antenna; then, the angular resolution of the lens can be expressed as sin θ m = mλ/D y , where λ, and θ m symbolize the carrier wavelength, and the angular position of the mth antenna relative to the x-axis, respectively.Let the central antenna be the reference element; then, the antennas can be indexed as −M, . . ., 0, . . ., M, where M is the half number of total elements.Furthermore, we also have θ m ≤ π/2; thus, we can assume that M = D y /λ without loss of generality [10].In this paper, we consider an LAA system composed of M = 2M + 1 array elements to receive farfield signals.According to [1, Lemma 1], the array response of the mth antenna to a signal from direction θ is given by a m (θ) = √ ψsinc (m − D y sin θ/λ), where ψ denotes the normalized effective aperture.
We assume L narrowband signals from K sources arriving at the LAA from distinct directions θ l , for l = 1, . . ., L. Let L k be the number of signals from the kth source, then we have of first and last signals departing from the kth source, respectively [4].Furthermore, a multipath signal s i (t), which is the replica of the direct signal s l k,1 (t) from the kth source, can be expressed as , where α i denotes the complex-valued attenuation coefficient.The received signal at the mth antenna at time t over T snapshots are given as where v m (t) represents the zero-mean Gaussian noise at the mth antenna.The entire received signal y (t) = [y −M (t) , . . .y 0 (t) , . .

. , y M (t)]
T ∈ C M ×1 can be rewritten in matrix form as y (t) = As (t) + v (t) where T ∈ C M ×1 denotes the noise vector.We assume that the noise and signals are mutually independent.In addition, the noise vector v (t) follows the circularly-symmetric Gaussian distribution CN 0, σ 2 v I M .Given the received signal vector, the signal covariance matrix (SCM) of y (t) is evaluated as R = E y (t) y H (t) .Similar to the results in [4, Eq. ( 5)], the (m, n)th element of R ∈ C M × M can be expressed as where δ m,n equals 1 if m = n, and 0 otherwise, and the pseudosignal s l,n can be expressed as where σ 2 s,k is the power of the direct signal from the kth transmitter, i.e., σ 3), it can be observed that the correlation of two signals transmitted from the kth source is given as s p,n / s q,n =α p /α q , for p = q ∈ [l k,1 , l k,2 ]; therefore, this yields the rank loss of R. In addition, the Vandermonde structures of the modified steering matrices that occur in the output of spatial smoothing algorithms [4], [9] cannot be formed by the array response of an antenna in an LAA since it follows a sinc function.As a result, the smoothed covariance matrices of the LAA have the full rank (M + 1) > L, making it impossible to discriminate between the signal and noise eigenvalues [4], [14], [15].Therefore, the high-resolution MUSIC algorithm cannot be applied for the LAA-aided multipath DoA estimation.

B. Conventional OMP Algorithm
In the COMP-based DoA estimation algorithm [11], a measurement matrix Ξ = a (θ 1 ) , a (θ 2 ) , . . ., a θ L ∈ R M × L contains L steering vectors.Let γ denote the resolution of DoA estimation in the angle range of [θ min , θ max ), then the set of DoA samples, and value of L can be expressed as Φ = {θ min , θ min + γ, . . ., θ max − γ}, and L = |Φ| = (θ max − θ min )/γ, respectively.The noisy received signals can then be rewritten as y (t) = Ξς (t) + v (t).The OMP can be employed to reconstruct the sparse signal ς (t) ∈ C L×1 by solving the optimization problem given as where ε represents the error-tolerant factor.Algorithm 1 describes the steps of the COMP scheme employed to solve the minimization problem (4).COMP incrementally constructs the support set ϑ by utilizing the matching step.This step involves a total of L − i + 1 2 M − 1 operations, entailing dot products between the steering vector a (θ l ) and Initialize i = 1, [0] = y (t), ϑ = ∅, and S = 0 L×T .

15:
end if end while 17: end for C M ×1 during the ith iteration.The subsequent aspect of Step 4, focused on identifying the largest absolute value, necessitates 2 L − i + 1 operations.Subsequently, the least square solution in Step 7 can be achieved by implementing the Cholesky or QR factorization methods [16].The total number of FLOPs in this step is approximately 2 Mi2 .Once the estimated signal ς [i] (t) is obtained, the residual is then updated as outlined in Step 9, which requires 2 Mi FLOPs, before being utilized to check the loop termination condition in Step 10.Let I t be the number of iterations needed to recover ς (t) = ς [It] (t).The entire complexity of the COMP to construct the matrix of recovered signals S = [ ς (1) , . . ., ς (T )] can be expressed as The COMP-sparse signal estimation is substantially degraded in a highly contaminated environment, as the noise elements become integrated into the reconstructed signal ς (t) [17].When the number of signals L is unknown, the termination condition based on ε cannot be executed effectively in a low-SNR environment [5].This leads to performance degradation and significant computational complexity.

A. Signal Covariance Matrix-Based OMP Method
Based on (2), the SCM R can be reformulated as In addition, it is worth highlighting that for any antenna index m and DoA θ satisfying we have a m (θ) ≈ √ ψ.This observation implies that the mth column of pseudosignal matrix S contains dominant values.On the other hand, other columns of S, whose indices m meet the condition |m − D y sin θ/λ| 0, have entries s l,m , which can be approximately zeros.It is noted that only a few entries in S have significantly larger amplitudes than others; thus, the mth column of R, denoted as R :,m , where m satisfies (7), is eligible for reconstructing the sparse where M is the set of active antenna indices, each satisfying (7).We note that the COMP algorithm can be reused to represent the SCM-OMP by modifying the input and output data.Specifically, the implementation of the SCM-OMP algorithm can be expressed as S = COMP R :,m m∈M , Ξ, ε , where COMP (•) represents procedure shown in Algorithm 1.
In addition, y (t), T , and S in COMP scheme are substituted by the selected column of noiseless SCM R :,m , the number of active antennas |M|, and the matrix of reconstructed pseudosignals S = ξ m m∈M ∈ C L×|M| , respectively.In the SCM-OMP algorithm, the SCM R is practically obtained by R = T t=1 x (t) x H (t) T .This step requires M 2 + M T FLOPs.In addition, the output spectrum of the SCM-OMP can be obtained by calculating Γ = diag S S H |M| .This   (7).However, from (2), we can observe that a signal arriving from a specific direction θ l , which is focused by the electromagnetic lens onto the mth antenna, satisfies (7) and results in an energy value of . On the other hand, the remaining diagonal elements of matrix R that correspond to inactive antennas are approximated as zeros.Consequently, the spectrum obtained by extracting the diagonal elements of R provides sharp peaks representing the indices of the antennas capturing the energy of incoming signals.In this scenario, some classical sharp peak recognition techniques, such as constant false alarm rate (CFAR) detector [18], or power spectrum-based convolutional neural network (PSCNet) [17], can be applied to identify M.
It is worth noting that the sharp peak recognition schemes [17], [18] only provide an index for the antenna that receives the most energy from a set of adjacent elements.This results in a misdetection in a case where two signals arrive at two consecutive antennas within a subarray, adversely affecting the performance of the OMP method.In addition, since condition ( 7

B. Multiple Sub-Covariance Matrices-Based OMP Algorithm
It can be observed in ( 9) that the computational complexity of the SCM-OMP scheme depends significantly on the number of active antennas.Specifically, when |M| is larger than T , the SCM-OMP exhibits higher complexity than the COMP method.However, it becomes apparent that the pseudosignal ξ m can be recovered using only a few entries from the mth column of R, given that R :,m is sparse.Based on Remark 1, it is obvious that the positions around the high-signal-energy diagonal entries of R :,m also exhibit significant magnitudes.Hence, we propose utilizing submatrices comprising some neighboring elements associated with dominant components on the spectrum d to facilitate pseudosignal reconstruction 2 .The cth submatrices can be then defined as Algorithm 2 summarizes the spatial spectrum reconstruction procedure based on the proposed MSCM-OMP scheme.It 2 We note that the set Mc may contain more than one index representing antennas with strong energy reception.Combining the two overlapping neighboring regions of the two peak indices prevents the repetition of applying the OMP algorithm to the intersection area.Thus, this reduces the complexity of the entire proposed algorithm.In addition to the case shown in Fig. 1(a), where two neighboring regions have no overlap, the two submatrices with intersecting elements can be merged into a larger submatrix, as visualized in Fig. 1(b).
is observed that Algorithm 2 constructs the entire output spectrum Γ based on C submatrices.Similar to (6), the computation of R c in Step 3 is given by where x c (t) where I c,m denotes the iteration number to recover ξ c,m , and It is clear that and ; therefore, based on ( 5), (9), and ( 12), our proposed MSCM-OMP algorithm requires a lower complexity compared to the aforementioned CS schemes.Subsequently, we consider the following lemma.These signals are transmitted from K = 3 sources and are recorded after T = 10 uniquely spaced time snapshots.In addition, we assume that the set of active antennas is correctly detected.In this section, we provide simulation results that offer comparative analyses between the aforementioned OMP-based DoA estimation schemes and other approaches, including the FTMR-supported MUSIC [4], ANM method [8], and sum signal covariance vector-based OMP (SSCV-OMP) [5].Furthermore, we utilize the root-mean-square-error (RMSE) metric to evaluate the performance of spectral-based algorithms.The RMSE metric with Q = 500 trials of the Monte Carlo simulation is defined as Fig. 2 illustrates the RMSE performance of our proposed algorithms in comparison with other spectral-based schemes.It can be observed that the performance gain of the MSCM-OMP method significantly increases with the increase in the number of antennas M or SNR.Fig. 2(a) demonstrates that the MSCM-OMP algorithm substantially outperforms other methods for M ≥ 30.It is worth highlighting that in Fig. 2(b), both the proposed SCM-and MSCM-OMP methods achieve better performance than other spectral-based DoA estimation methods in a contaminated environment, especially at SNR ≤ 5 dB.For example, the RMSE of SCM-and MSCM-OMP algorithms are 4 • and 0.6 • at SNR = 0 dB, respectively.In contrast, the error of the high-resolution MUSIC scheme is 10 • , which is approximately 2.5 and 17 times higher than those of the SCM-OMP and MSCM-OMP, respectively.In addition, Fig. 2(c) shows that the MSCM-OMP achieves the best performance gain with any values of T .Moreover, the error of the MSCM-OMP decreases as the number of snapshots T increases.
Fig. 3 provides a comparison of the computational complexity among several DoA estimation methods, including the FTMR-supported MUSIC, ANM scheme, and OMP-based DoA estimators such as COMP, SSCV-OMP, SCM-OMP, and MSCM-OMP algorithms.It is shown that the number of computations of the MSCM-OMP scheme is lower than that of other methods for M ≥ 30.In addition, the complexity of the MSCM-OMP decreases significantly as M increases.
This work was supported by the Australian Research Council Discovery Project under Grant DP220101158.(Corresponding author: Trong-DaiHoang.) where the element in the (l, n) position of the pseudosignal matrix S ∈ C L× M is provided in (3).It is evident that the element r m,n , where m = n, remains uncontaminated, given the mutual independence of signals and noise.Noise components are strictly concentrated along the diagonal of matrix R, and σ 2 v can be interpreted as the minimum value among the entries {r m,m } M m=−M .Next, we define r min Δ = min {r m,m } M m=−M .Then, the noiseless covariance matrix R Δ = AS can be approximated as

Fig. 1 :
Fig. 1: Formations of submatrices when neighbors of two detected peak indices (a) have no common elements, and (b) share an intersection.

Δ=
x mc,1:m c,|Mc | (t) ∈ C |Mc|×1 .This step requires |M c | 2 + |M c | T operations.The total complexity of the MSCM-OMP primarily hinges on the iterations of constructing S c = ξ c,m m∈Mc in Step 6.Furthermore, the cth region of the entire spatial spectrum Γ can be separately obtained based on S c , as shown in Step 7. Finally, the total complexity of the MSCM-OMP can be expressed as

Lemma 1 .
In a cluster of active antennas, the coverage range of candidate DoAs becomes narrower as the size of LAA or the number of antennas increases.Proof.The DoA coverage function of the cth antenna cluster with respect to the normalized lens dimension D y /λ can be defined as f Dy λ Δ = sin −1 m c,|Mc | Dy/λ − sin −1 mc,1 Dy/λ .The first-order derivative of f (D y /λ) can be calculated as df (Dy/λ) d(Dy/λ) = 1 Dy/λ g (m c,1 ) − g m c,|Mc| , where the function g (x) is defined as g (x) Δ = x √ (Dy/λ) 2 −x 2 .It is worth highlighting that g (x) increases monotonically since dg(x) dx > 0 for any value of x.Therefore, we have df (Dy/λ) d(Dy/λ) < 0. In other words, the function f (D y /λ) decreases monotonically as the D y -λ ratio increases.Moreover, given M = D y /λ , the number of DoA grid points is reduced as M or D y /λ becomes larger.This concludes the proof.By consideringLemma 1,and (12), it is worth emphasizing that the complexity of MSCM-OMP decreases significantly as the number of antennas increases.The operation number of the MSCM-OMP converges to OC c=1 |M c | 2 |Φ c | + |Φ c | 2 |M c |as M approaches infinity.Specifically, both |M c |, and |Φ c | decrease substantially as M becomes larger.In contrast, the spatial smoothing-aided MUSIC[4], the ANM method[8], COMP, SSCV-OMP[5], and SCM-OMP algorithms require the com-plexity of O M 4 , O M 4 , O M LT , O M 2 T , andO M L |M| + M 2 T as Δ m → ∞,respectively.This implies that the MSCM-OMP scheme is more suitable for estimating signal directions with a large-scale LAA than other existing methods.IV.SIMULATION RESULTS AND ANALYSISFor model settings, the LAA systems with an effective aperture of ψ = 1 and half-array size of M ∈[10, 100]   are utilized to receive signals from distinct directions in a range of θ ∈ [−60 • , 60 • ).The resolution of DoA samples is set to γ = 0.1 • , resulting in an output spatial spectrum length of L = 1200.For the subsequent simulations, we focus on the LAA with M = 30, which receives L = 5 signals from directions θ = {−37.4• , 17.6 • , −50.8 • , −19.8 • , 42.2 • }.
and ε.Output: S ∈ C M ×T .