# Compressed Sensing PARALIND Decomposition-Based Coherent Angle Estimation for Uniform Rectangular Array

Authors: Wu Wei; Xu Le; Zhang Xiaofei; Li Jianfeng
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

In this paper, the topic of coherent two-dimensional direction of arrival (2D-DOA) estimation is investigated. Our study jointly utilizes the compressed sensing (CS) technique and the parallel profiles with linear dependencies (PARALIND) model and presents a 2D-DOA estimation algorithm for coherent sources with the uniform rectangular array. Compared to the traditional PARALIND decomposition, the proposed algorithm owns lower computational complexity and smaller data storage capacity due to the process of compression. Besides, the proposed algorithm can obtain autopaired azimuth angles and elevation angles and can achieve the same estimation performance as the traditional PARALIND, which outperforms some familiar algorithms presented for coherent sources such as the forward backward spatial smoothing-estimating signal parameters via rotational invariance techniques (FBSS-ESPRIT) and forward backward spatial smoothing-propagator method (FBSS-PM). Extensive simulations are provided to validate the effectiveness of the proposed CS-PARALIND algorithm.


---

## Body

## 1. Introduction

Array signal processing has aroused considerable concerns in recent decades owing to its extensive engineering application in satellite communication, radar, sonar, and some other fields [1–4]. In array signal processing, spectrum estimation, also known as direction of arrival (DOA) estimation, is a crucial issue. Till now, there are already many neoteric algorithms [5–8] proposed for DOA estimation with linear array, which include estimating signal parameters via rotational invariance techniques (ESPRIT) algorithm [5, 6], multiple signal classification (MUSIC) algorithm [7], and propagator method (PM) [8]. Compared with linear array, rectangular array can measure both azimuth angle and elevation angle, and hence 2D-DOA estimation with rectangular array has motivated enormous investigations. Many traditional DOA estimation methods have already been extended to 2D-DOA estimation [9–12]. For example, the 2D-ESPRIT algorithms [9–11] have utilized the invariance property to obtain the 2D-DOA estimation with uniform rectangular array, and the 2D-MUSIC [12] algorithm has also proved to be applicable. Some other algorithms such as angle estimation with generalized coprime planar array [13] and 2D-DOA estimation with nested subarrays [14] have better performance for 2D-DOA estimation as well.However, in many practical situations, the complex propagation environment usually leads to the presence of coherent signals, which makes it complicated to obtain the valid DOA estimation. Therefore, the research on coherent angle estimation has gained great significance. The algorithms mentioned above are only applicable to noncoherent sources and the coherent sources will lead to severe invalidity for these methods. Some other coherent estimation methods including traditional forward spatial smoothing (FSS) or forward backward spatial smoothing (FBSS) [15, 16] algorithm have good estimation performance. Exploiting coprime multiple-input multiple-output (MIMO) radar, [17] has proposed the DOA estimation method for mixed coherent and uncorrelated targets and [18] obtained the DOA estimation of coherent sources via fourth-order cumulants.The trilinear decomposition, namely, parallel factor (PARAFAC) technique [19–21], has been widely employed to resolve the problem of angular estimation with rectangular array [22]. However, the decomposition of PARAFAC model fails to work if it contains coherent sources. The parallel profiles with linear dependencies (PARALIND) model [23, 24] discussed in this paper can be regarded as a generalization of PARAFAC model and is efficient to solve the problem of coherent DOA estimation, where the PARAFAC method usually cannot present significant results. In [25, 26], the PARALIND decomposition method has already been successfully applied to obtain the coherent sources estimation with MIMO radar and acoustic vector-sensor array. However, the traditional algorithms based on PARALIND decomposition involve heavy computational burden and huge capacity consumption of data storage.Compressed sensing (CS) [27, 28] has aroused considerable concern, which is introduced to areas including channel estimation, image, beamforming, and radar [29–32]. Specifically, the angular information of sources can be structured as a sparse vector and hence the CS technique can be directly utilized [33, 34]. By applying CS theory, many novel parameter estimation algorithms have been proposed for different scenarios. Reference [35] proposed a CS-based angle estimation method for noncircular sources which can be applied to MIMO radar, while [33, 34] combined the CS technique with PARAFAC model and proposed an estimation algorithm for joint angle and frequency. Based on the CS PARAFAC model, a 2D-DOA estimation method was presented in [36] with a uniform rectangular array, whereas it works only for noncoherent sources.In this paper, we propose a CS-PARALIND algorithm for coherent sources to extract the DOA estimates with a uniform rectangular array by generalizing the work presented in [36] and the CS theory for PARALIND decomposition. We first construct the received data model which can be transformed to the PARALIND model. Then, we compress the data model and perform the PARALIND decomposition on it to achieve the estimation of compressed direction matrices. Finally, to acquire the 2D-DOA estimation, we formulate a sparse recovery problem which can be solved by the orthogonal matching pursuit (OMP) method [37]. Note that the compression process only compresses the directional matrix and the source matrix while the coherent matrix remains the same. Therefore, the coherent structure of impinging signals is not destroyed after compression. Comparing to [36], we summarize the main contributions of our research in this paper:( 1 ) We construct the received data model of coherent signals for uniform rectangular array which is suitable for the PARALIND decomposition.( 2 ) We generalize the method in [36] and propose the compressed sensing PARALIND (CS-PARALIND) model. Notably, there is no existing report for CS-PARALIND decomposition so far as we know.( 3 ) We develop the CS-PARALIND decomposition method to obtain the coherent 2D-DOA estimation for uniform rectangular array and give the detailed description of it.In this paper, the proposed algorithm can obtain autopaired 2D-DOA estimation of coherent signals. In addition, the corresponding correlated matrix can also be obtained. The proposed algorithm can attain the same angle estimation performance as the traditional PARALIND method [25]. Compared to the FBSS-ESPRIT and the FBSS-PM algorithm, our method can achieve better DOA estimation performance. In addition, due to the process of compression, the proposed algorithm consumes lower computational burden and requires limited storage capacity in practical application. The Cramér-Rao bound (CRB) for the DOA estimation with uniform rectangular array is also provided in this paper. A series of simulation results verify the effectiveness of our approach.The remainder of our paper is organized as follows: Section2 presents the received data model of coherent signals with uniform rectangular array. Section 3 depicts the detailed derivation of the proposed CS-PARALIND algorithm as well as the uniqueness certification and complexity analysis. Numerical simulations are exhibited in Section 4, and we conclude this paper in Section 5.Notation. ⊕, ∘, and ⊗, respectively, representHadamard product,Khatri–Rao product, andKronecker product. (.)H, (.)T, (.)∗, (.)-1, and (.)+ stand for the operations of conjugate-transpose, transpose, conjugation, inverse, and pseudoinverse. TheFrobenius norm andl0–norm are denoted by ∙F and ∙0. diag(b) is a diagonal matrix composed of elements in vector b while diag-1(A) produces a column vector composed of the diagonal elements of matrix A. Dn(A) denotes a diagonal matrix which consists of then-th row of A. span(H) represents the subspace spanned by the columns of H. vec(A) means stacking the columns of matrix A. ∂ means partial derivatives.


## 2. Data Model

Assume that there areK far-field narrow-band signals with DOA (θk,φk) impinging on a uniform rectangular array consisting of M×N sensors, where φk is the azimuth angle of the k-th signal, θk represents the elevation angle, and the distance between any two adjacent elements is d. The sources numberK is known and the noise is additive white Gaussian which is assumed to be independent and uncorrelated with received signals. Figure 1 shows the structure of the array.Figure 1
The structure of uniform rectangular array.For the first subarray in the uniform rectangular array, the received data att-th time can be represented as [36] (1)x1t=Axs0t+n1t,where Ax=[ax(θ1,φ1),ax(θ2,φ2),⋯,ax(θK,φK)] denotes the directional matrix of the first subarray with ax(θk,φk)=[1,exp(j2πdcos⁡φksin⁡θk/λ),⋯,exp(j2πM-1dcos⁡φksin⁡θk/λ)]T and j=-1. The received noise of the first subarray is denoted by n1(t), and s0(t) is the K×1 source vector of theK signals. Then it follows that the received data of n-th subarray at t-th time can be written as [36](2)xnt=AxΦn-1s0t+nnt,where Φ=diag(ej2πdsin⁡θ1sin⁡φ1/λ,⋯,ej2πdsin⁡θKsin⁡φK/λ) and nn(t) stands for the received noise of the n-th subarray. Therefore, the received data of the entire rectangular array can be denoted by [38](3)xt=x1tx2t⋮xNt=AxAxΦ⋮AxΦN-1s0t+n1tn2t⋮nNt,or more compactly as(4)xt=Ay∘Axs0t+nt,where Ay=[ay(θ1,φ1),ay(θ2,φ2),⋯,ay(θK,φK)] with ay(θk,φk)=[1,exp(j2πdsin⁡φksin⁡θk/λ),⋯,exp(j2πN-1dsin⁡φksin⁡θk/λ)] and n(t) is the noise of the whole rectangular array with n(t)=[n1(t)T,n2(t)T,…,nN(t)T]T.By exploitingJ snapshots, the received data can be denoted as [36](5)X~=xt1,xt2,⋯,xtJ,where X~ stands for the noisy received signals, and rewriting (5) in matrix form, we obtain [38](6)X~=X~1X~2⋮X~N=AxD1AyAxD2Ay⋮AxDNAyS0+N1N2⋮NN=Ay∘AxS0+N,where S0=[s0(t1),s0(t2),⋯,s0(tJ)]∈CK×J and X~n=AxDn(Ay)S0+Nn. The noise ofJ snapshots is N=[n(t1),n(t2),⋯,n(tJ)]∈CMN×J.The received data matrix shown in (6) is a traditional PARAFAC model [19–22], the decomposition of which is usually nonunique if there exist coherent signals [19]. Assume that among theKreceived sources, there are K1 groups of coherent signals; then (6) can be represented as [25](7)X~=Ay∘AxΓS+N,where S∈CK1×J is the source matrix of J snapshots of K1 noncoherent signals. Γ∈CK×K1 is the corresponding correlated matrix satisfying ΓS=S0.


## 3. CS-PARALIND Decomposition-Based Algorithm

The PARALIND algorithm suffers from expensive computational cost, especially in the case of large number of sensors or snapshots. Specifically, the angular information of received sources can be structured as a sparse vector and hence the CS technique can be directly utilized [33]. To counter this problem, we bring in the CS theory by compressing the received signal in (7) to a smaller matrix, followed by the PARALIND decomposition and, finally, acquire the 2D-DOA estimates by exploiting sparse recovery method.


### 3.1. Compression

We compress the received dataX_∈CM×J×N into a smaller matrix X′_∈CM′×J′×N′with three compression matrices U∈CM×M′, V∈CN×N′, and W∈CJ×J′, with M′<M, J′<J and N′<N, which is illustrated in Figure 2.Figure 2
Compression processing.The three compression matrices can be constructed via some random special matrices,e.g., the random Gaussian, Bernoulli, and partial Fourier matrices or the Tucker3 decomposition introduced in [33, 39].ApplyingU, V, and W to (7), we can obtain the compressed received data X~′∈C(M′N′×J′) as [27](8)X~′=VT⊗UTX~MN×JW=VT⊗UTAy∘AxΓSW+VT⊗UTNW.And according to the property of Khatri–Rao product [33], we have (VT⊗UT)(Ay∘Ax)=(VTAy)∘(UTAx), and then (8) can be further simplified to(9)X~′=Ay′∘Ax′ΓS′+N′,where Ax′(M′×K)=UTAx, Ax′(M′×K)=UTAx, Ay′(N′×K)=VTAy, S′(K×J′)=SW, and N′=(VT⊗UT)NW.


### 3.2. PARALIND Decomposition

To transform (9) to PARALIND model, we define Y~=X~′T, or equivalently, (10)Y~=ΓS′TAy′∘Ax′T+N′T=Y+N′T,where Y=(ΓS′)T(Ay′∘Ax′)T denotes the noise-free compressed received matrix.Obviously, (10) is equivalent to the PARALIND model [20] and the least fitting for it amounts to [23](11)minΓ,S,Ay′,Ax′Y~-ΓS′TAy′∘Ax′TF2.Under noise-free condition, (10) can be rewritten as(12)Y=Y1,Y2,⋯,YN′=ΓS′TD1Ay′Ax′T,ΓS′TD2Ay′Ax′T,⋯,ΓS′TDN′Ay′Ax′T,where Yn=(ΓS′)TDnAy′Ax′T∈CJ′×M′, n=1,2,⋯,N′.After vectoringYn, we have [38](13)vecYn=vecΓS′TDnAy′Ax′T=Ax′DnAy′⊗S′TvecΓT.Stacking these vectors by columns(14)vecY1vecY2⋮vecYN′=Ax′D1Ay′⊗S′TAx′D2Ay′⊗S′T⋮Ax′DN′Ay′⊗S′TvecΓT.Equation (14) also can be denoted as (15)vecY=Ay′∘Ax′⊗S′TvecΓT.Then under noise environment, vec(ΓT) can be obtained by(16)vecΓT=Ay′∘Ax′⊗S′T+vecY~.Now we can update the value of Γ via transforming vec(ΓT) to its original modality.According to (11), the least square (LS) update for S′T is (17)S′T=Y~Ax′∘Ay′ΓT+.From (12), it follows that the covariance matrix of Yn can be written as(18)∑n=1N′YnHYn=∑n=1N′YnHΓS′TDnAy′Ax′T=Ax′∗∑n=1N′Dn∗Ay′ΓS′∗ΓS′TDnAy′Ax′T.As Ax′ is of full column rank and ∑n=1N′Dn∗(Ay′)(ΓS′)∗(ΓS′)TDn(Ay′) is nonsingular [24], we can obtain Ax′∗ via(19)Ax′∗=∑n=1N′Y~nHΓS′TDnAy′∑n=1N′Dn∗Ay′ΓS′∗ΓS′TDnAy′-1=∑n=1N′Y~nHΓS′TDnAy′∑n=1N′ΓS′∗ΓS′T⊕Ay′HAy′-1.SubstituteYn of (12) into the following (20)ΓS′∗YnAx′=ΓS′∗ΓS′TDnAy′Ax′TAx′.Taking the diagonal elements of both sides in (20), we have(21)diag-1ΓS′∗YnAx′=ΓS′∗ΓS′T⊕Ax′TAx′diag-1DnAy′.By enforcing left multiplication operation of (21), we get(22)diag-1DnAy′=ΓS′∗ΓS′T⊕Ax′TAx′-1diag-1ΓS′∗Y~nAx′.And hence, Ay can be easily obtained form (22).For the PARALIND decomposition above, according to (17), (18), and (19), we repeatedly update the estimation of Ax′, Ay′, S′, and Γ until convergence. Define the sum of squared residuals (SSR) of thek-th iteration as SSRk=∑i=1N′M′∑j=1J′cij2,where cij is the (i,j) element of C=Y~-(Γ^S^′)T(A^y′∘A^x′)T. Define the convergence rate of PARALIND decomposition as SSRrate=(SSRk-SSRk-1)/SSRk-1 [40]. When SSRrate is smaller than a certain small value (determined by the noise level), then the above iteration process can be considered convergent.Based on the PARALIND decomposition, the estimates ofAx′ and Ay′ can be achieved by(23)A^x′=Ax′ΠΔx+Wx=UTAxΠΔx+Wx,(24)A^y′=Ay′ΠΔy+Wy=VTAyΠΔy+Wy,where Π stands for the permutation matrix and Δy and Δx are the diagonal scaling matrix of Ay′ and Ax′. In addition, Wx and Wy denote the estimation error.Remark 1.
In the proposed algorithm, the scale ambiguity can be eliminated by direct normalization while the permutation ambiguity makes no difference in the angle estimation.


### 3.3. DOA Estimation with Sparsity

By utilizingA^x′, A^y′, the 2D-DOA estimation can be obtained with sparsity. With noiseless case, we use a^xk′ and a^yk′ to denote the k-th column of A^x′, A^y′, respectively. From (23) and (24), it follows that [36](25)a^xk′=UHρxkaxk,(26)a^yk′=VHρykayk,where axk and ayk represent the k-th column in the directional matrices Ax and Ay. In addition, ρxk and ρyk denote the scaling coefficients.Then we construct two Vandermonde matricesAsx∈CM×Q and Asy∈CN×Q (Q>>M,Q>>N), the columns of which consist of the steering vectors corresponding to each potential source location [36].(27)Asx=asx1,asx2,⋯,asxQ=11⋯1ej2πdg1/λej2πdg2/λ⋮ej2πdgQ/λ⋮⋮⋮⋮ej2πM-1dg1/λej2πM-1dg2/λ⋯ej2πM-1dgQ/λ,(28)Asy=asy1,asy2,⋯,asyQ=11⋯1ej2πdg1/λej2πdg2/λ⋮ej2πdgQ/λ⋮⋮⋮⋮ej2πN-1dg1/λej2πN-1dg2/λ⋯ej2πN-1dgQ/λ,where g is a sampling vector with g(q)=-1+2q/Q, q=0,1,…,Q. Asx and Asy can be referred to as the overcomplete dictionary for our 2D-DOA estimation [36]. Then (25), (26) can be converted to(29)a^xk′=UTAsxxsk,k=1,⋯,K,(30)a^yk′=VTAsyysk,k=1,⋯,K,where xsk and ysk are sparse, the estimates of which can be cast as an optimization problem, subject to l0-norm constraint.(31)mina^xk′-UHAsxxskF2,st.xsk0=1,(32)mina^yk′-VHAsyyskF2,st.ysk0=1.xsk and ysk can be gained by the OMP recovery method [37].Denote the maximum modulus of elements inxsk and ysk as qxk and qyk, respectively. Then the corresponding elements g(qxk) and g(qyk) in Asx and Asy are just the sin⁡θkcos⁡φk and sin⁡θksin⁡φk estimation. Define rk=g(qxk)+jg(qyk) and the estimates of elevation angles and azimuth angles can be obtained by (33)θ^k=sin-1absrk,k=1,2,…,K,(34)φ^k=anglerk,k=1,2,…,K,where abs(.) outputs the modulus value and angle(.) computes the angular part of a complex number. Finally, the autopaired estimates of elevation and azimuth angles can be attained by employing the paired relationship in the directional matrices obtained by PARALIND decomposition.


### 3.4. The Procedure of the Proposed Algorithm

Till now, the angle estimation of received coherent signals with a uniform rectangular array has been acquired and we provide the major steps as follows:Step 1.
Compress the received dataX~, construct the PARALIND model Y~ according to (8), and then initialize the value of Ax′, Ay′, Γ, and S′.Step 2.
According to (16), (17), (19), and (22), repeatedly update the estimates of Ax′, Ay′, S′, and Γ according to the convergence conditions.Step 3.
According to (31)-(34), the elevation and azimuth angles can be estimated by Asx and Asy.


## 3.1. Compression

We compress the received dataX_∈CM×J×N into a smaller matrix X′_∈CM′×J′×N′with three compression matrices U∈CM×M′, V∈CN×N′, and W∈CJ×J′, with M′<M, J′<J and N′<N, which is illustrated in Figure 2.Figure 2
Compression processing.The three compression matrices can be constructed via some random special matrices,e.g., the random Gaussian, Bernoulli, and partial Fourier matrices or the Tucker3 decomposition introduced in [33, 39].ApplyingU, V, and W to (7), we can obtain the compressed received data X~′∈C(M′N′×J′) as [27](8)X~′=VT⊗UTX~MN×JW=VT⊗UTAy∘AxΓSW+VT⊗UTNW.And according to the property of Khatri–Rao product [33], we have (VT⊗UT)(Ay∘Ax)=(VTAy)∘(UTAx), and then (8) can be further simplified to(9)X~′=Ay′∘Ax′ΓS′+N′,where Ax′(M′×K)=UTAx, Ax′(M′×K)=UTAx, Ay′(N′×K)=VTAy, S′(K×J′)=SW, and N′=(VT⊗UT)NW.


## 3.2. PARALIND Decomposition

To transform (9) to PARALIND model, we define Y~=X~′T, or equivalently, (10)Y~=ΓS′TAy′∘Ax′T+N′T=Y+N′T,where Y=(ΓS′)T(Ay′∘Ax′)T denotes the noise-free compressed received matrix.Obviously, (10) is equivalent to the PARALIND model [20] and the least fitting for it amounts to [23](11)minΓ,S,Ay′,Ax′Y~-ΓS′TAy′∘Ax′TF2.Under noise-free condition, (10) can be rewritten as(12)Y=Y1,Y2,⋯,YN′=ΓS′TD1Ay′Ax′T,ΓS′TD2Ay′Ax′T,⋯,ΓS′TDN′Ay′Ax′T,where Yn=(ΓS′)TDnAy′Ax′T∈CJ′×M′, n=1,2,⋯,N′.After vectoringYn, we have [38](13)vecYn=vecΓS′TDnAy′Ax′T=Ax′DnAy′⊗S′TvecΓT.Stacking these vectors by columns(14)vecY1vecY2⋮vecYN′=Ax′D1Ay′⊗S′TAx′D2Ay′⊗S′T⋮Ax′DN′Ay′⊗S′TvecΓT.Equation (14) also can be denoted as (15)vecY=Ay′∘Ax′⊗S′TvecΓT.Then under noise environment, vec(ΓT) can be obtained by(16)vecΓT=Ay′∘Ax′⊗S′T+vecY~.Now we can update the value of Γ via transforming vec(ΓT) to its original modality.According to (11), the least square (LS) update for S′T is (17)S′T=Y~Ax′∘Ay′ΓT+.From (12), it follows that the covariance matrix of Yn can be written as(18)∑n=1N′YnHYn=∑n=1N′YnHΓS′TDnAy′Ax′T=Ax′∗∑n=1N′Dn∗Ay′ΓS′∗ΓS′TDnAy′Ax′T.As Ax′ is of full column rank and ∑n=1N′Dn∗(Ay′)(ΓS′)∗(ΓS′)TDn(Ay′) is nonsingular [24], we can obtain Ax′∗ via(19)Ax′∗=∑n=1N′Y~nHΓS′TDnAy′∑n=1N′Dn∗Ay′ΓS′∗ΓS′TDnAy′-1=∑n=1N′Y~nHΓS′TDnAy′∑n=1N′ΓS′∗ΓS′T⊕Ay′HAy′-1.SubstituteYn of (12) into the following (20)ΓS′∗YnAx′=ΓS′∗ΓS′TDnAy′Ax′TAx′.Taking the diagonal elements of both sides in (20), we have(21)diag-1ΓS′∗YnAx′=ΓS′∗ΓS′T⊕Ax′TAx′diag-1DnAy′.By enforcing left multiplication operation of (21), we get(22)diag-1DnAy′=ΓS′∗ΓS′T⊕Ax′TAx′-1diag-1ΓS′∗Y~nAx′.And hence, Ay can be easily obtained form (22).For the PARALIND decomposition above, according to (17), (18), and (19), we repeatedly update the estimation of Ax′, Ay′, S′, and Γ until convergence. Define the sum of squared residuals (SSR) of thek-th iteration as SSRk=∑i=1N′M′∑j=1J′cij2,where cij is the (i,j) element of C=Y~-(Γ^S^′)T(A^y′∘A^x′)T. Define the convergence rate of PARALIND decomposition as SSRrate=(SSRk-SSRk-1)/SSRk-1 [40]. When SSRrate is smaller than a certain small value (determined by the noise level), then the above iteration process can be considered convergent.Based on the PARALIND decomposition, the estimates ofAx′ and Ay′ can be achieved by(23)A^x′=Ax′ΠΔx+Wx=UTAxΠΔx+Wx,(24)A^y′=Ay′ΠΔy+Wy=VTAyΠΔy+Wy,where Π stands for the permutation matrix and Δy and Δx are the diagonal scaling matrix of Ay′ and Ax′. In addition, Wx and Wy denote the estimation error.Remark 1.
In the proposed algorithm, the scale ambiguity can be eliminated by direct normalization while the permutation ambiguity makes no difference in the angle estimation.


## 3.3. DOA Estimation with Sparsity

By utilizingA^x′, A^y′, the 2D-DOA estimation can be obtained with sparsity. With noiseless case, we use a^xk′ and a^yk′ to denote the k-th column of A^x′, A^y′, respectively. From (23) and (24), it follows that [36](25)a^xk′=UHρxkaxk,(26)a^yk′=VHρykayk,where axk and ayk represent the k-th column in the directional matrices Ax and Ay. In addition, ρxk and ρyk denote the scaling coefficients.Then we construct two Vandermonde matricesAsx∈CM×Q and Asy∈CN×Q (Q>>M,Q>>N), the columns of which consist of the steering vectors corresponding to each potential source location [36].(27)Asx=asx1,asx2,⋯,asxQ=11⋯1ej2πdg1/λej2πdg2/λ⋮ej2πdgQ/λ⋮⋮⋮⋮ej2πM-1dg1/λej2πM-1dg2/λ⋯ej2πM-1dgQ/λ,(28)Asy=asy1,asy2,⋯,asyQ=11⋯1ej2πdg1/λej2πdg2/λ⋮ej2πdgQ/λ⋮⋮⋮⋮ej2πN-1dg1/λej2πN-1dg2/λ⋯ej2πN-1dgQ/λ,where g is a sampling vector with g(q)=-1+2q/Q, q=0,1,…,Q. Asx and Asy can be referred to as the overcomplete dictionary for our 2D-DOA estimation [36]. Then (25), (26) can be converted to(29)a^xk′=UTAsxxsk,k=1,⋯,K,(30)a^yk′=VTAsyysk,k=1,⋯,K,where xsk and ysk are sparse, the estimates of which can be cast as an optimization problem, subject to l0-norm constraint.(31)mina^xk′-UHAsxxskF2,st.xsk0=1,(32)mina^yk′-VHAsyyskF2,st.ysk0=1.xsk and ysk can be gained by the OMP recovery method [37].Denote the maximum modulus of elements inxsk and ysk as qxk and qyk, respectively. Then the corresponding elements g(qxk) and g(qyk) in Asx and Asy are just the sin⁡θkcos⁡φk and sin⁡θksin⁡φk estimation. Define rk=g(qxk)+jg(qyk) and the estimates of elevation angles and azimuth angles can be obtained by (33)θ^k=sin-1absrk,k=1,2,…,K,(34)φ^k=anglerk,k=1,2,…,K,where abs(.) outputs the modulus value and angle(.) computes the angular part of a complex number. Finally, the autopaired estimates of elevation and azimuth angles can be attained by employing the paired relationship in the directional matrices obtained by PARALIND decomposition.


## 3.4. The Procedure of the Proposed Algorithm

Till now, the angle estimation of received coherent signals with a uniform rectangular array has been acquired and we provide the major steps as follows:Step 1.
Compress the received dataX~, construct the PARALIND model Y~ according to (8), and then initialize the value of Ax′, Ay′, Γ, and S′.Step 2.
According to (16), (17), (19), and (22), repeatedly update the estimates of Ax′, Ay′, S′, and Γ according to the convergence conditions.Step 3.
According to (31)-(34), the elevation and azimuth angles can be estimated by Asx and Asy.


## 4. Performance Analysis

### 4.1. Complexity Analysis of the Proposed Algorithm

We analyze the complexity of the proposed algorithm in this subsection. For the proposed algorithm, the computational complexity isO(M′N′MNJ+JM′N′J′+n1(M′N′J′2K2K12+2KK1+2K+K1 + M′N′2K12+K2+KK1+2K + K2J′N′+J′+N′2+N′+M′+2+K3K13 + 2K3+K13+KK1J′N′+1))+KPM+N, where n1 denotes the time of iterations. With regard to the traditional PARALIND decomposition [23], it has (n2(MNJ(2K2K12+2KK1+2K+K1)+MN(2K12+K2+KK1+2K) + K2(JN+J+N2+N+M+2) + K3K13+2K3+K13+KK1J(N+1)) + 2K2(M+N)+6K2), where n2 represents the time of iterations. In the comparison of computational complexity, we have the parameters of M=N=10, K=3, and K1=2. In addition, assume that n1=n2=20, P=200, and M′/M=N′/N=J′/J=0.5. The comparison of complexity versus J is depicted in Figure 3, which concludes that the computational burden of the proposed algorithm is decreased compared to the traditional PARALIND algorithm.Figure 3
Comparison of complexity versus different snapshots (J).


### 4.2. CRB

DefineA=[ay(θ1,φ1)⊗ax(θ1,φ1),⋯,ay(θK,φK)⊗ax(θK,φK)] and according to [41], we derive the CRB(35)CRB=σ22JReDHΠA⊥D⊕P^sT-1,where J is the number of snapshots and ak is the k-th column of A. σ2 denotes the noise power. ΠA⊥=IMN-AAHA-1AH and IMN is an MN×MN identity matrix. D= [∂a1/∂θ1,∂a2/∂θ2,...,∂aK/∂θK,∂a1/∂φ1,∂a2/∂φ2,...,∂aK/∂φK] and P^s=1/J∑t=1JstsHt.


### 4.3. Advantages

We summarize the advantages of the proposed algorithm as follows:( 1 ) The proposed algorithm possesses lower computational cost and requires smaller data storage, due to the compression operation.( 2 ) The proposed algorithm can be applied to the coherent signals. Furthermore, the corresponding correlated matrix also can be obtained.( 3 ) The proposed algorithm can achieve the same angle estimation performance with the traditional PARALIND method and outperforms the FBSS-ESPRIT and FBSS-PM algorithm.


## 4.1. Complexity Analysis of the Proposed Algorithm

We analyze the complexity of the proposed algorithm in this subsection. For the proposed algorithm, the computational complexity isO(M′N′MNJ+JM′N′J′+n1(M′N′J′2K2K12+2KK1+2K+K1 + M′N′2K12+K2+KK1+2K + K2J′N′+J′+N′2+N′+M′+2+K3K13 + 2K3+K13+KK1J′N′+1))+KPM+N, where n1 denotes the time of iterations. With regard to the traditional PARALIND decomposition [23], it has (n2(MNJ(2K2K12+2KK1+2K+K1)+MN(2K12+K2+KK1+2K) + K2(JN+J+N2+N+M+2) + K3K13+2K3+K13+KK1J(N+1)) + 2K2(M+N)+6K2), where n2 represents the time of iterations. In the comparison of computational complexity, we have the parameters of M=N=10, K=3, and K1=2. In addition, assume that n1=n2=20, P=200, and M′/M=N′/N=J′/J=0.5. The comparison of complexity versus J is depicted in Figure 3, which concludes that the computational burden of the proposed algorithm is decreased compared to the traditional PARALIND algorithm.Figure 3
Comparison of complexity versus different snapshots (J).


## 4.2. CRB

DefineA=[ay(θ1,φ1)⊗ax(θ1,φ1),⋯,ay(θK,φK)⊗ax(θK,φK)] and according to [41], we derive the CRB(35)CRB=σ22JReDHΠA⊥D⊕P^sT-1,where J is the number of snapshots and ak is the k-th column of A. σ2 denotes the noise power. ΠA⊥=IMN-AAHA-1AH and IMN is an MN×MN identity matrix. D= [∂a1/∂θ1,∂a2/∂θ2,...,∂aK/∂θK,∂a1/∂φ1,∂a2/∂φ2,...,∂aK/∂φK] and P^s=1/J∑t=1JstsHt.


## 4.3. Advantages

We summarize the advantages of the proposed algorithm as follows:( 1 ) The proposed algorithm possesses lower computational cost and requires smaller data storage, due to the compression operation.( 2 ) The proposed algorithm can be applied to the coherent signals. Furthermore, the corresponding correlated matrix also can be obtained.( 3 ) The proposed algorithm can achieve the same angle estimation performance with the traditional PARALIND method and outperforms the FBSS-ESPRIT and FBSS-PM algorithm.


## 5. Simulation Results

Suppose thatK=3 signals with two coherent signals and one noncoherent signal impinge on a rectangular array, where the correlated matrix is Γ=101010T. In the following simulations, we exploit the root mean square error (RMSE) to evaluate the DOA estimation performance and it is presented by (36)RMSE=1K∑k=1K1L∑l=1Lφ^k,l-φk2+θ^k,l-θk2,where φ^k,l, θ^k,l are the estimations of φk, θk in l-th simulation, and the times of Monte-Carlo simulations are indicated by L. In the following simulations, we set L=1000.The locations of the three sources are(θ1,φ1)=(15∘,10∘), (θ2,φ2)=(25∘,30∘), and (θ3,φ3)=(35∘,50∘). M,N,J,K are the number of rows and columns of the rectangular array, snapshots, and sources, respectively. And assume that the rectangular array is uniform and the distance of any two adjacent sensors is λ/2.Simulation 1.
Figure4 shows the 2D-DOA estimation results of our proposed algorithm with SNR=0dB and SNR=20dB, respectively. The size of PARALIND model is defined as M×N×J. The received data matrix has the size of 10×10×200 in this simulation, which becomes 7×7×140(M′×N′×J′) after compression. Figure 4 indicates that the proposed algorithm is effective for 2D-DOA estimation of coherent sources.Figure 4
DOA estimation of the proposed algorithm with (a) SNR=0dB, (b) SNR=20dB.
(a) (b)Simulation 2.
In Figure5, we exhibit the comparison of the DOA estimation performance among the proposed algorithm, FSS-ESPRIT, FSS-PM, FBSS-ESPRIT, and FBSS-PM algorithm, and the traditional PARALIND algorithm. The received data matrix has the size of 10×10×100 and is compressed to the size of 8×8×80. It is illustrated clearly in Figure 5 that the approach algorithm can achieve the same DOA estimation performance compared with the traditional PARALIND algorithm and, furthermore, outperforms the other four algorithms.Figure 5
DOA estimation performance comparison.Simulation 3.
The DOA estimation performance result of the proposed algorithm versus compression ratio is provided in Figure6, where the compression ratio is P=M′/M=N′/N=J′/J. In this part, the size of the original PARALIND model is 15×15×200 with P=0.4, P=0.6, and P=0.8. It is shown in Figure 6 that the estimation performance of our method improves withPgetting larger.Figure 6
DOA estimation performance under different compression ratio.Simulation 4.
The angle estimation performance of the proposed algorithm versusJ is illustrated in Figure 7, where the rectangular array is structured as M=10, N=10 and M′=N′=8, J′/J=0.8. It is verified in Figure 7 that the proposed algorithm can obtain improved DOA estimation accuracy with a larger snapshots.Figure 7
DOA estimation performance under different snapshots.Simulation 5.
Figure8 shows DOA estimation performance of the proposed algorithm with varied M. Assume that N=N′=10 and J=J′=200 with M′/M=2/3. It is attested by Figure 8 that the estimation performance of the proposed method gets better with the increasing M.Figure 8
DOA estimation performance under differentM.Simulation 6.
The estimation performance of the proposed method versusN is illustrated in Figure 9, where M=15, J=200 with M′/M=J′/J=N′/N=0.8. It is indicated obviously that the proposed method obtains improved performance with larger N.Figure 9
DOA estimation performance under differentN.


## 6. Conclusions

We jointly utilize the CS theory and the PARALIND decomposition in this paper and propose a CS-PARALIND algorithm with a uniform rectangular array to extract the DOA estimates of coherent signals. The proposed algorithm can attain autopaired 2D-DOA estimates, and benefiting from the compression procedure, lower computational cost, and smaller demand for storage capacity can be achieved. Extensive simulations corroborate that the proposed approach obtains superior estimation performance compared to the traditional PARALIND algorithm and, significantly, outperforms the FBSS-PM and FBSS- ESPRIT algorithm.



---

*Source: 3619818-2019-06-24.xml*