# Real-Valued Weighted Subspace Fitting Algorithm for DOA Estimation with Block Sparse Recovery

Authors: Liangliang Li; Xianpeng Wang; Jinmei Shi; Xiang Lan
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

In this paper, the problem of direction-of-arrival (DOA) estimation for strictly noncircular sources under the condition of unknown mutual coupling is concerned, and then a robust real-valued weighted subspace fitting (WSF) algorithm is proposed via block sparse recovery. Inspired by noncircularity, the real-valued coupled extended array output with double array aperture is first structured via exploiting the real-valued conversion. Then, an efficient real-valued block extended sparse recovery model is constructed by performing the parameterized decoupling operation to avoid the unknown mutual coupling and noncircular phase effects. Thereafter, the WSF framework is investigated to recover the real-valued block sparse matrix, where the spectrum of real-valued NC MUSIC-like is utilized to design a weighted matrix for strengthening the solutions sparsity. Eventually, DOA estimation is achieved based on the support set of the reconstructed block sparse matrix. Owing to the combination of noncircularity, parametrized decoupling thought, and reweighted strategy, the proposed method not only effectively achieves high-precision estimation, but also efficiently reduces the computational complexity. Plenty of simulation results demonstrate the effectiveness and efficiency of the proposed method.


---

## Body

## 1. Introduction

Thanks to the growing maturity of array signal processing technology, parameter estimation gradually occupies an important position in the fields of vehicle positioning, radar, medical diagnosis, and so on [1]. As one of the bases for parameter estimation, direction-of-arrival (DOA) estimation has been a hot topic for decades accompanied by a series of work [2–4]. Afterwards, benefiting from the increasing development of multiple-input multiple-output (MIMO) technique, MIMO radar architectures have been developed to provide high degrees of freedom (DOF) and resolution for DOA estimation [5]. Unfortunately, the distance between sensors decreases as the number of antennas increases for a fixed array aperture. It means that it is quite possible for closely-spaced sensors to suffer from the unknown mutual coupling effect. Thereby, it is worthwhile to study DOA estimation with strong robustness. In this way, this paper mainly investigates the robust DOA estimation of strictly noncircular sources with unknown mutual coupling.Generally speaking, many DOA estimation attempts can be roughly divided into subspace-based methods [6–9] and sparse signal recovery (SSR) methods [10–13]. Multiple signal classification (MUSIC) method [6] uses the decoupled noise subspace to first achieve super-resolution direction finding, different from estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [7] based on the decoupled signal subspace. It should be pointed out that these approaches have difficulty in achieving satisfactory performance under low signal-to-noise ratio (SNR), insufficient snapshots, or correlated sources. Thereafter, sparse recovery technique offers a feasible perspective to overcome these drawbacks, which can be categorized into norm optimization estimators [10, 11] and sparse Bayesian learning (SBL) approaches [12, 13]. Furthermore, it has been demonstrated that SSR algorithms are better than subspace-based methods in challenging circumstances, such as unsatisfactory SNR or inadequate snapshots [14].It can be found that the above methods study circular signals by default. However, in recent years, DOA estimation of noncircular sources has received extensive attention in parameter estimation [15]. This is largely due to its wide distribution and natural superiorities. To the best of our knowledge, noncircular sources are commonly seen in practical communication systems [16], such as amplitude modulation (AM) and binary phase shift keying (BPSK). More importantly, noncircular sources can achieve higher accuracy and detect more targets than the default circular sources [17]. Subsequently, numerous algorithms [18–24] have been presented for noncircular sources that show the advantage in accuracy. On the one hand, there are lots of attempts achieved by subspace technology [18–20]. As shown in [18], noncircular MUSIC (NC MUSIC) algorithm is derived via combining noncircularity with MUSIC principle. Whereas, large-scale spectral peak search results in relatively high computational complexity. After that, noncircular root MUSIC (NC Root-MUSIC) approach [19] and noncircular conjugate ESPRIT (NC C-ESPRIT) algorithm [20] are introduced for tackling the above problem. On the other hand, DOA estimation for noncircular sources is implemented from the perspective of sparse reconstruction [21–24]. In [21, 22], the joint sparsity-aware schemes for array and monostatic MIMO radar system are put forward, respectively. With the in-depth research on sparsity, not only block sparsity but also rank sparsity are simultaneously utilized to model a nuclear norm penalty framework for enhancing the solutions sparsity [23]. Although this method has superiorities in estimation accuracy and resolution, it is computationally expensive. Thereby, a unitary nuclear norm minimization strategy [24] is further presented to reduce the computational complexity.It is noted that the array manifolds of the above methods are normally assumed to be ideal. Nevertheless, such hypothesis may not be applicable to practice due to the existence of array manifold perturbations, like mutual coupling [25, 26]. It is generally believed that there may be unknown mutual coupling between closely-spaced antennas affected by the interaction of space electromagnetic fields [26]. This perturbation leads to undesired array manifold, thereby degrading or even invalidating the estimation performance of these approaches. Afterwards, a large number of calibration ideas are designed to deal with the problem of unknown mutual coupling [27–38].For one thing, a series of calibrations [27–34] for circular sources have been attempted to estimate DOAs. In [27], the unknown mutual coupling is modeled as a complex band symmetric Toeplitz structure, and then additional auxiliary sensors are added to compensate. Similarly, the selection matrix is further designed by setting the antennas at both ends of the original array to be auxiliary sensors [28]. Unfortunately, these approaches can only maintain normal direction finding at the expense of array aperture. For preserving the array aperture as much as possible, the parameterized decoupling idea [29] is introduced to decouple the angle parameter and mutual coupling coefficients. Although this method uses whole data, its application scope is still limited because it belongs to subspace-based methods. Different from these efforts using subspace technology, relevant works [30–34] on sparse recovery have also been carried out. As introduced in [30], a revised l1-SVD (singular value decomposition) algorithm is structured by designing a specific selection matrix in array. Analogously, the selection matrix is further implanted into the MIMO framework [31]. Actually, both the auxiliary sensors and the selection matrix can be regarded as two embodiments of array compensation. But they all sacrifice the array aperture. Aiming at this drawback, an effective block sparse recovery (BSR) approach [32] is presented by replacing array compensation with parameterized decoupling. Moreover, a reweighted BSR algorithm [33] and a weighted subspace fitting (WSF) method [34] are further reported for acquiring higher accuracy.For another, some studies [35–40] on noncircular sources have been done to estimate DOAs. In [35], a selection matrix is first constructed to remove the negative influence so as to directly apply ESPRIT principle. Similar to [27, 28], it is achieved at the expense of array aperture. Subsequently, an efficient real-valued rank reduction method [36] using MUSIC principle is derived to effectively avoid the unknown mutual coupling effect and protect the precious array aperture. However, these methods are still subject to the limitations of subspace technology, unlike robust SSR algorithms. In view of this shortcoming, the joint reweighted sparsity-inducing scheme based on SVD principle [37] and WSF principle [38] are put forward, respectively. Whereas, their computational complexity is relatively higher than that of subspace-based methods in [35, 36].In this work, an efficient real-valued WSF algorithm with block sparse recovery is presented for DOA estimation of strictly noncircular sources under unknown mutual coupling. First, a real-valued block extended sparse recovery model is formed to avoid the unknown mutual coupling and noncircular phase effects. Subsequently, the regularization framework between sparsity penalty and subspace fitting error is investigated. Finally, a real-valued reweighted block sparse recovery approach is explored to achieve WSF for DOA estimation. The proposed method effectively maintains high accuracy and efficiently reduces the computational load. The simulation results confirm the correctness of the above deduce.The main contributions of the proposed method are summarized as follows:(a)
Perform a real-valued conversion to reduce the computational burden, and then construct a real-valued coupled extended data by exploiting noncircularity.(b)
Eliminate the unknown mutual coupling and noncircular phase interferences through parameterized decoupling operation without array aperture loss.(c)
Structure a real-valued noncircular MUSIC-like (NC MUSIC-like) weighted matrix to enhance the solutions sparsity.(d)
Develop a robust real-valued WSF framework to estimate DOAs by block sparse recovery.It is noted that some important notations adopted in this article are defined in Table1.Table 1
Some important notations.
NotationsDefinitions⋅T , ⋅∗ and ⋅HTranspose, conjugate, and conjugate-transpose operations⋅† and ⋅Pseudo-inverse and absolute value operationsRe⋅ and Im⋅Real and imaginary part operations0M×KM×K dimensional zero matrixIMM×M dimensional identity matrixdiag⋅ and blkdiag⋅Diagonalization and block diagonalization operationsE⋅Mathematical expectation operationtr⋅ and det⋅Trace and determinant operations⋅0 , ⋅1, and ⋅2l0-norm, l1-norm, and l2-norm⋅FFrobenius norm


## 2. Data Model for DOA Estimation

### 2.1. Problem Formulation

Suppose thatK far-field uncorrelated narrowband NC sources sd,kk=1K incident on a uniform linear array (ULA) equipped with M omnidirectional antennas. The distinct DOAs can be denoted as θ=θ1,θ2,…,θK. Then, the array output in the ideal environment can be structured as(1)xt=Asdt+nt,where xt=x1t,x2t,…,xMtT∈ℂM×1 is the received data. nt=n1t,n2t,…,nMtT∈ℂM×1 stands for the complex additive Gaussian white noise vector with zero mean. Meanwhile, sdt=sd,1t,sd,2t,…,sd,KtT∈ℂK×1 means the noncircular signal vector. A=aθ1,aθ2,…,aθK∈ℂM×K indicates the ideal array manifold matrix. As the k th column of matrix A, aθk is known as array manifold corresponding to k th target and satisfies aθk=1,ρθk,…,ρM−1θkT∈ℂM×1, where ρθk=ejνθk with νθk=−2πd/λnsinθk. d represents the distance between adjacent antennas and λn denotes the signal wavelength. Obviously, the data model in (1) is not affected by any array manifold perturbations, such as mutual coupling [37] and gain-phase error [9]. Unfortunately, as the number of antennas increases, the distance between sensors decreases. Hence, due to the fact that space electromagnetic field interacts with each other, closely-spaced antennas are vulnerable to unknown mutual coupling, as illustrated in Figure 1.Figure 1
Mutual coupling model of ULA.It has been demonstrated that the mutual coupling coefficients between antennas are inversely proportional to their spacing. That is to say, the larger the sensors spacing, the weaker the mutual coupling effect, and the smaller the corresponding coefficients. Moreover, when the antennas are far enough away from each other, it is reasonable to ignore the impact of unknown mutual coupling. Then, a structure of complex banded symmetric Toeplitz in [27] is utilized to model the mutual coupling matrix (MCM), i.e.,(2)G=Toeplitz1,g1,…,gH−1,01×M−H=1g1…gH−1g11g1…gH−10⋮⋱⋱⋱⋱⋱gH−1⋯g11g1⋯gH−1⋱⋮⋱⋱⋱⋱⋱gH−1⋯g11g1⋯gH−10⋱⋮⋱⋱⋱⋮gH−1⋯g11g1gH−1⋯g11,where gcc=1H−1 refer to the unknown nonzero mutual coupling coefficients, whose c th element gc=λcejφc is made up of amplitude coefficient λc and phase coefficient φc, respectively. It can be clearly seen that there are H nonzero mutual coupling coefficients, which satisfy 0<gH−1<gH−2<,…,<g2<g1<g0=1. Then, the ideal array manifold is affected by the unknown mutual coupling in antennas, which should be revised as(3)a∧θk=Gaθk.Thereby, the practical array output under the condition of unknown mutual coupling can be written as(4)yt=GAsdt+nt,where yt=y1t,y2t,…,yMtT∈ℂM×1 represents the actual received data disturbed by unknown mutual coupling, unlike xt in (1).According to the above introduction,sdt denotes the noncircular source, which means that its noncircular rate ξ ranges from 0 to 1, including the upper limit [38]. When referring to the maximum noncircular rate ξ=1, the radiation signal can be defined as the strictly noncircular source, like AM modulation signal. In this paper, strictly noncircular sources are considered. It has revealed in [22] that the complex strictly NC source sdt in (4) can be further expressed as(5)sdt=Φst,where Φ=diagejϕ1,ejϕ2,…,ejϕK∈ℂK×K stands for the rotation phase shift matrix corresponding to ϕ=ϕ1,ϕ2,…,ϕK, which can be arbitrary for each source. st=s1t,s2t,…,sKtT∈ℝK×1 is the real-valued signal vector corresponding to the complex-valued vector sdt. In this way, taking (5) back to (4), the actual array output can be represented as(6)yt=GAsdt+nt=GAΦst+nt.


### 2.2. Real-Valued Conversion for Noncircular Sources

In view of the noncircularity advantages, many researches on noncircular sources directly construct the extended signal model achieved by the received data and its conjugate form. However, the data belongs to the complex domain, which inevitably leads to the high computational burden. Different from the above classical processing in the complex domain, a real-valued conversion is first applied to the received data for structuring an extended data model in the real domain [36]. Thanks to the real-valued transformation, the computation load is greatly reduced to accelerate the direction finding speed. Then, following the idea in [36], the real and imaginary parts of the actual array output can be extracted as(7)yRt=Reyt=yt+y∗t2=GAΦ+G∗A∗Φ∗2st+Rent=ARst+nRt,(8)yIt=Imyt=yt−y∗t2j=GAΦ−G∗A∗Φ∗2jst+Imnt=AIst+nIt,where nRt∈ℝM×1 and nIt∈ℝM×1 express the real and imaginary components achieved by the complex noise vector nt∈ℂM×1. AR=aRθ1,ϕ1,G,aRθ2,ϕ2,G,…,aRθK,ϕK,G∈ℝM×K and AI=aIθ1,ϕ1,G,aIθ2,ϕ2,G,…,aIθK,ϕK,G∈ℝM×K are the virtual coupled array manifold matrices of yRt and yIt, respectively. Moreover, as one of the columns in AR and AI, aRθ,ϕ,G and aIθ,ϕ,G simultaneously contain the unknown mutual coupling coefficients and noncircular phase. They can be represented as(9)aRθ,ϕ,G=Gaθejϕ+G∗A∗θe−jϕ2,aIθ,ϕ,G=Gaθejϕ−G∗A∗θe−jϕ2j.Combining (7) and (8), a real-valued coupled extended signal model can be designed as(10)yRIt=yRtyIt​=​ARAIst+nRtnIt=ARIst+nRIt,where ARI=aRIθ1,ϕ1,G,aRIθ2,ϕ2,G,…,aRIθK,ϕK,G∈ℝ2M×K is the real-valued coupled extended steering matrix. Each column in ARI denotes the coupled extended array manifold and takes the following structure:(11)aRIθ,ϕ,G=aRθ,ϕ,GaIθ,ϕ,G.Then, the covariance matrix ofyRIt can be written as(12)R∧=EyRItyRIHt=ARIRsARIH+σ2I2M,where σ2 expresses the corresponding noise power. Rs=EstsHt indicates the signal covariance matrix, whose rank K′ rests with the source correlation. This paper assumes K′=K because of the uncorrelated sources presupposition. In fact, R∧ can only be obtained when the number of snapshots approaches infinity. However, it is clearly unavailable and eventually replaced by its maximum likelihood estimation R in reality. R can be computed by finite snapshots T, which takes the following form:(13)R=1T∑t=1TyRItyRIHt,where R refers to the sampling covariance matrix. Then, applying eigenvalue decomposition to R, yields(14)R=∑m=12MλmδmδmH=EsΩsEsH+EnΩnEnH,where λmm=12M mean the eigenvalues and satisfy λ1≥λ2≥,…,≥λK>λK+1=,…,=λ2M. δmm=12M are the eigenvectors corresponding to the eigenvalues λmm=12M. What is more, K larger eigenvalues and their corresponding eigenvectors are utilized to formulate the diagonal matrix Ωs and the signal subspace Es, respectively, i.e., Ωs=diagλ1,λ2,…,λK∈ℝK×K and Es=δ1,δ2,…,δK∈ℝ2M×K. In like manner, the diagonal matrix Ωn=diagλK+1,λK+2,…,λ2M∈ℝ2M−K×2M−K and the noise subspace En=δK+1,δK+2,…,δ2M∈ℝ2M×2M−K are composed of 2M−K smaller eigenvalues and the corresponding eigenvectors [38].As introduced in [39], the steering matrix spans the same range subspace as the signal subspace. Similarly, the signal subspace Es lies in the range space of the real-valued coupled extended steering matrix ARI, which indicates that Es and ARI satisfy(15)Es=ARIU.where U denotes a column full rank matrix with K×K dimension. Unfortunately, it is hard for (15) to estimate U due to the coexistence of unknown parameters, such as mutual coupling coefficients and noncircular phase. Thereby, a robust estimator should be designed to overcome these disturbances.


### 2.3. Parameterized Decoupling Operation

It is noted that the real-valued coupled extended array manifold in (11) is affected by unknown mutual coupling and noncircular phase, resulting in the failure of many existing ideal direction finding algorithms. Inspired by [36], the parameterized decoupling thought in the real domain is exploited to deal with the above problem.Through parameterizing the virtual coupled array manifoldsaRθ,ϕ,G and aIθ,ϕ,G, yields(16)aRθ,ϕ,G=Gaθejϕ+G∗a∗θe−jϕ2=Ψ⌢θ,ϕ,GT⌢θΣ⌢θ,ϕ,G−Ψ⌣θ,ϕ,GT⌣θΣ⌣θ,ϕ,G,(17)aIθ,ϕ,G=Gaθejϕ−G∗a∗θe−jϕ2j=Ψ⌣θ,ϕ,GT⌢θΣ⌣θ,ϕ,G+Ψ⌢θ,ϕ,GT⌣θΣ⌢θ,ϕ,G,where(18)T⌢θ=1cosνθ⋱0cosH−1νθ⋮0cosM−Hνθ⋱cosM−1νθ,T⌣θ=0sinνθ⋱0sinH−1νθ⋮0sinM−Hνθ⋱sinM−1νθ,Ψ⌢θ,ϕ,G=∑c=1−HH−1λccosφc+ϕ+c⋅νθ,Ψ⌣θ,ϕ,G=∑c=1−HH−1λcsinφc+ϕ+c⋅νθ,Σ⌢θ,ϕ,G=τ⌢1θ,…,τ⌢H−1θ,1,σ⌢1θ,…,σ⌢H−1θT,Σ⌣θ,ϕ,G=τ⌣1θ,…,τ⌣H−1θ,1,σ⌣1θ,…,σ⌣H−1θT,τ⌢hθ=Ψ⌢θ,ϕ,G−∑c=hH−1λccosφc+ϕ−c⋅νθΨ⌢θ,ϕ,G,τ⌣hθ=Ψ⌣θ,ϕ,G−∑c=hH−1λcsinφc+ϕ−c⋅νθΨ⌣θ,ϕ,G,σ⌢hθ=Ψ⌢θ,ϕ,G−∑c=H−hH−1λccosφc+ϕ+c⋅νθΨ⌢θ,ϕ,G,σ⌣hθ=Ψ⌣θ,ϕ,G−∑c=H−hH−1λcsinφc+ϕ+c⋅νθΨ⌣θ,ϕ,G,where h=1,2,…,H−1 . T⌢θ∈ℝM×2H−1 and T⌣θ∈ℝM×2H−1 are the real-valued block matrices and only depend on angle information. For briefness, F=2H−1 is defined in what follows.It is worth emphasizing thatΨ⌢θ,ϕ,G and Ψ⌣θ,ϕ,G stand for two constants. They rely on three parameters, i.e., angle parameter, mutual coupling coefficients, and noncircular phase. Additionally, Ψ⌢θ,ϕ,G≠0 and Ψ⌣θ,ϕ,G≠0 occur with extremely high probability except for a few very special circumstances. Thereby, this article defaults Ψ⌢θ,ϕ,G≠0 and Ψ⌣θ,ϕ,G≠0.Then, bringing (16) and (17) back to (11), aRIθ,ϕ,G∈ℝ2M×1 can be decoupled as(19)aRIθ,φ,G=aRθ,φ,GaIθ,φ,G=Ψ⌢θ,φ,GT⌢θΣ⌢θ,φ,G−Ψ⌣θ,φ,GT⌣θΣ⌣θ,φ,GΨ⌣θ,φ,GT⌢θΣ⌣θ,φ,G+Ψ⌢θ,φ,GT⌣θΣ⌢θ,φ,G=T⌢θ−T⌣θT⌣θT⌢θ⏟T¯θΨ⌢θ,φ,GΣ⌢θ,φ,GΨ⌣θ,φ,GΣ⌣θ,φ,G⏟Λθ,φ,G,where T¯θ∈ℝ2M×2F and Λθ,ϕ,G∈ℝ2F×1 are block matrix and block vector, respectively. It can be discovered that the real-valued coupled extended array manifold aRIθ,ϕ,G∈ℝ2M×1 is decoupled into two parts: T¯θ and Λθ,ϕ,G. T¯θ only rests with angle parameter. Therefore, it can be seemed as a new decoupled extended steering vector, similar to aRIθ,ϕ,G∈ℝ2M×1 in (11). While Λθ,ϕ,G subjects to the interferences of unknown mutual coupling and noncircular phase.According to (19), the coupled signal model in (10) can be further decoupled as (20)yRIt=yRtyIt=ARAIst+nRtnIt=ARIst+nRIt=A∧RIΔst+nRIt=A∧RIsRIt+nRIt,where(21)A^RI=T¯θ1,T¯θ2,…,T¯θK∈R2M×2FK,Δ=blkdiagΛθ1,φ1,G,Λθ2,φ2,G,…,ΛθK,φK,G,where A∧RI denotes the real-valued block extended array manifold matrix formed by T¯θkk=1,2,…,K. It separates angle parameter from disturbance factors, such as mutual coupling coefficients and noncircular phase, making it only depend on DOAs information. The block diagonal matrix Δ∈ℝ2FK×K is combined with the real signal vector st to construct a novel block signal vector sRIt. i.e., sRIt=Δst∈ℝ2FK×1. It can be found that the 2Fk−2F+1 th to 2Fk th rows in sRIt correspond to k th element in st. Besides, both the new extended array manifold matrix A∧RI and the corresponding signal vector sRIt in (20) have block structures for each target, unlike the original coupled extended steering matrix ARI and the signal vector st in (10).


## 2.1. Problem Formulation

Suppose thatK far-field uncorrelated narrowband NC sources sd,kk=1K incident on a uniform linear array (ULA) equipped with M omnidirectional antennas. The distinct DOAs can be denoted as θ=θ1,θ2,…,θK. Then, the array output in the ideal environment can be structured as(1)xt=Asdt+nt,where xt=x1t,x2t,…,xMtT∈ℂM×1 is the received data. nt=n1t,n2t,…,nMtT∈ℂM×1 stands for the complex additive Gaussian white noise vector with zero mean. Meanwhile, sdt=sd,1t,sd,2t,…,sd,KtT∈ℂK×1 means the noncircular signal vector. A=aθ1,aθ2,…,aθK∈ℂM×K indicates the ideal array manifold matrix. As the k th column of matrix A, aθk is known as array manifold corresponding to k th target and satisfies aθk=1,ρθk,…,ρM−1θkT∈ℂM×1, where ρθk=ejνθk with νθk=−2πd/λnsinθk. d represents the distance between adjacent antennas and λn denotes the signal wavelength. Obviously, the data model in (1) is not affected by any array manifold perturbations, such as mutual coupling [37] and gain-phase error [9]. Unfortunately, as the number of antennas increases, the distance between sensors decreases. Hence, due to the fact that space electromagnetic field interacts with each other, closely-spaced antennas are vulnerable to unknown mutual coupling, as illustrated in Figure 1.Figure 1
Mutual coupling model of ULA.It has been demonstrated that the mutual coupling coefficients between antennas are inversely proportional to their spacing. That is to say, the larger the sensors spacing, the weaker the mutual coupling effect, and the smaller the corresponding coefficients. Moreover, when the antennas are far enough away from each other, it is reasonable to ignore the impact of unknown mutual coupling. Then, a structure of complex banded symmetric Toeplitz in [27] is utilized to model the mutual coupling matrix (MCM), i.e.,(2)G=Toeplitz1,g1,…,gH−1,01×M−H=1g1…gH−1g11g1…gH−10⋮⋱⋱⋱⋱⋱gH−1⋯g11g1⋯gH−1⋱⋮⋱⋱⋱⋱⋱gH−1⋯g11g1⋯gH−10⋱⋮⋱⋱⋱⋮gH−1⋯g11g1gH−1⋯g11,where gcc=1H−1 refer to the unknown nonzero mutual coupling coefficients, whose c th element gc=λcejφc is made up of amplitude coefficient λc and phase coefficient φc, respectively. It can be clearly seen that there are H nonzero mutual coupling coefficients, which satisfy 0<gH−1<gH−2<,…,<g2<g1<g0=1. Then, the ideal array manifold is affected by the unknown mutual coupling in antennas, which should be revised as(3)a∧θk=Gaθk.Thereby, the practical array output under the condition of unknown mutual coupling can be written as(4)yt=GAsdt+nt,where yt=y1t,y2t,…,yMtT∈ℂM×1 represents the actual received data disturbed by unknown mutual coupling, unlike xt in (1).According to the above introduction,sdt denotes the noncircular source, which means that its noncircular rate ξ ranges from 0 to 1, including the upper limit [38]. When referring to the maximum noncircular rate ξ=1, the radiation signal can be defined as the strictly noncircular source, like AM modulation signal. In this paper, strictly noncircular sources are considered. It has revealed in [22] that the complex strictly NC source sdt in (4) can be further expressed as(5)sdt=Φst,where Φ=diagejϕ1,ejϕ2,…,ejϕK∈ℂK×K stands for the rotation phase shift matrix corresponding to ϕ=ϕ1,ϕ2,…,ϕK, which can be arbitrary for each source. st=s1t,s2t,…,sKtT∈ℝK×1 is the real-valued signal vector corresponding to the complex-valued vector sdt. In this way, taking (5) back to (4), the actual array output can be represented as(6)yt=GAsdt+nt=GAΦst+nt.


## 2.2. Real-Valued Conversion for Noncircular Sources

In view of the noncircularity advantages, many researches on noncircular sources directly construct the extended signal model achieved by the received data and its conjugate form. However, the data belongs to the complex domain, which inevitably leads to the high computational burden. Different from the above classical processing in the complex domain, a real-valued conversion is first applied to the received data for structuring an extended data model in the real domain [36]. Thanks to the real-valued transformation, the computation load is greatly reduced to accelerate the direction finding speed. Then, following the idea in [36], the real and imaginary parts of the actual array output can be extracted as(7)yRt=Reyt=yt+y∗t2=GAΦ+G∗A∗Φ∗2st+Rent=ARst+nRt,(8)yIt=Imyt=yt−y∗t2j=GAΦ−G∗A∗Φ∗2jst+Imnt=AIst+nIt,where nRt∈ℝM×1 and nIt∈ℝM×1 express the real and imaginary components achieved by the complex noise vector nt∈ℂM×1. AR=aRθ1,ϕ1,G,aRθ2,ϕ2,G,…,aRθK,ϕK,G∈ℝM×K and AI=aIθ1,ϕ1,G,aIθ2,ϕ2,G,…,aIθK,ϕK,G∈ℝM×K are the virtual coupled array manifold matrices of yRt and yIt, respectively. Moreover, as one of the columns in AR and AI, aRθ,ϕ,G and aIθ,ϕ,G simultaneously contain the unknown mutual coupling coefficients and noncircular phase. They can be represented as(9)aRθ,ϕ,G=Gaθejϕ+G∗A∗θe−jϕ2,aIθ,ϕ,G=Gaθejϕ−G∗A∗θe−jϕ2j.Combining (7) and (8), a real-valued coupled extended signal model can be designed as(10)yRIt=yRtyIt​=​ARAIst+nRtnIt=ARIst+nRIt,where ARI=aRIθ1,ϕ1,G,aRIθ2,ϕ2,G,…,aRIθK,ϕK,G∈ℝ2M×K is the real-valued coupled extended steering matrix. Each column in ARI denotes the coupled extended array manifold and takes the following structure:(11)aRIθ,ϕ,G=aRθ,ϕ,GaIθ,ϕ,G.Then, the covariance matrix ofyRIt can be written as(12)R∧=EyRItyRIHt=ARIRsARIH+σ2I2M,where σ2 expresses the corresponding noise power. Rs=EstsHt indicates the signal covariance matrix, whose rank K′ rests with the source correlation. This paper assumes K′=K because of the uncorrelated sources presupposition. In fact, R∧ can only be obtained when the number of snapshots approaches infinity. However, it is clearly unavailable and eventually replaced by its maximum likelihood estimation R in reality. R can be computed by finite snapshots T, which takes the following form:(13)R=1T∑t=1TyRItyRIHt,where R refers to the sampling covariance matrix. Then, applying eigenvalue decomposition to R, yields(14)R=∑m=12MλmδmδmH=EsΩsEsH+EnΩnEnH,where λmm=12M mean the eigenvalues and satisfy λ1≥λ2≥,…,≥λK>λK+1=,…,=λ2M. δmm=12M are the eigenvectors corresponding to the eigenvalues λmm=12M. What is more, K larger eigenvalues and their corresponding eigenvectors are utilized to formulate the diagonal matrix Ωs and the signal subspace Es, respectively, i.e., Ωs=diagλ1,λ2,…,λK∈ℝK×K and Es=δ1,δ2,…,δK∈ℝ2M×K. In like manner, the diagonal matrix Ωn=diagλK+1,λK+2,…,λ2M∈ℝ2M−K×2M−K and the noise subspace En=δK+1,δK+2,…,δ2M∈ℝ2M×2M−K are composed of 2M−K smaller eigenvalues and the corresponding eigenvectors [38].As introduced in [39], the steering matrix spans the same range subspace as the signal subspace. Similarly, the signal subspace Es lies in the range space of the real-valued coupled extended steering matrix ARI, which indicates that Es and ARI satisfy(15)Es=ARIU.where U denotes a column full rank matrix with K×K dimension. Unfortunately, it is hard for (15) to estimate U due to the coexistence of unknown parameters, such as mutual coupling coefficients and noncircular phase. Thereby, a robust estimator should be designed to overcome these disturbances.


## 2.3. Parameterized Decoupling Operation

It is noted that the real-valued coupled extended array manifold in (11) is affected by unknown mutual coupling and noncircular phase, resulting in the failure of many existing ideal direction finding algorithms. Inspired by [36], the parameterized decoupling thought in the real domain is exploited to deal with the above problem.Through parameterizing the virtual coupled array manifoldsaRθ,ϕ,G and aIθ,ϕ,G, yields(16)aRθ,ϕ,G=Gaθejϕ+G∗a∗θe−jϕ2=Ψ⌢θ,ϕ,GT⌢θΣ⌢θ,ϕ,G−Ψ⌣θ,ϕ,GT⌣θΣ⌣θ,ϕ,G,(17)aIθ,ϕ,G=Gaθejϕ−G∗a∗θe−jϕ2j=Ψ⌣θ,ϕ,GT⌢θΣ⌣θ,ϕ,G+Ψ⌢θ,ϕ,GT⌣θΣ⌢θ,ϕ,G,where(18)T⌢θ=1cosνθ⋱0cosH−1νθ⋮0cosM−Hνθ⋱cosM−1νθ,T⌣θ=0sinνθ⋱0sinH−1νθ⋮0sinM−Hνθ⋱sinM−1νθ,Ψ⌢θ,ϕ,G=∑c=1−HH−1λccosφc+ϕ+c⋅νθ,Ψ⌣θ,ϕ,G=∑c=1−HH−1λcsinφc+ϕ+c⋅νθ,Σ⌢θ,ϕ,G=τ⌢1θ,…,τ⌢H−1θ,1,σ⌢1θ,…,σ⌢H−1θT,Σ⌣θ,ϕ,G=τ⌣1θ,…,τ⌣H−1θ,1,σ⌣1θ,…,σ⌣H−1θT,τ⌢hθ=Ψ⌢θ,ϕ,G−∑c=hH−1λccosφc+ϕ−c⋅νθΨ⌢θ,ϕ,G,τ⌣hθ=Ψ⌣θ,ϕ,G−∑c=hH−1λcsinφc+ϕ−c⋅νθΨ⌣θ,ϕ,G,σ⌢hθ=Ψ⌢θ,ϕ,G−∑c=H−hH−1λccosφc+ϕ+c⋅νθΨ⌢θ,ϕ,G,σ⌣hθ=Ψ⌣θ,ϕ,G−∑c=H−hH−1λcsinφc+ϕ+c⋅νθΨ⌣θ,ϕ,G,where h=1,2,…,H−1 . T⌢θ∈ℝM×2H−1 and T⌣θ∈ℝM×2H−1 are the real-valued block matrices and only depend on angle information. For briefness, F=2H−1 is defined in what follows.It is worth emphasizing thatΨ⌢θ,ϕ,G and Ψ⌣θ,ϕ,G stand for two constants. They rely on three parameters, i.e., angle parameter, mutual coupling coefficients, and noncircular phase. Additionally, Ψ⌢θ,ϕ,G≠0 and Ψ⌣θ,ϕ,G≠0 occur with extremely high probability except for a few very special circumstances. Thereby, this article defaults Ψ⌢θ,ϕ,G≠0 and Ψ⌣θ,ϕ,G≠0.Then, bringing (16) and (17) back to (11), aRIθ,ϕ,G∈ℝ2M×1 can be decoupled as(19)aRIθ,φ,G=aRθ,φ,GaIθ,φ,G=Ψ⌢θ,φ,GT⌢θΣ⌢θ,φ,G−Ψ⌣θ,φ,GT⌣θΣ⌣θ,φ,GΨ⌣θ,φ,GT⌢θΣ⌣θ,φ,G+Ψ⌢θ,φ,GT⌣θΣ⌢θ,φ,G=T⌢θ−T⌣θT⌣θT⌢θ⏟T¯θΨ⌢θ,φ,GΣ⌢θ,φ,GΨ⌣θ,φ,GΣ⌣θ,φ,G⏟Λθ,φ,G,where T¯θ∈ℝ2M×2F and Λθ,ϕ,G∈ℝ2F×1 are block matrix and block vector, respectively. It can be discovered that the real-valued coupled extended array manifold aRIθ,ϕ,G∈ℝ2M×1 is decoupled into two parts: T¯θ and Λθ,ϕ,G. T¯θ only rests with angle parameter. Therefore, it can be seemed as a new decoupled extended steering vector, similar to aRIθ,ϕ,G∈ℝ2M×1 in (11). While Λθ,ϕ,G subjects to the interferences of unknown mutual coupling and noncircular phase.According to (19), the coupled signal model in (10) can be further decoupled as (20)yRIt=yRtyIt=ARAIst+nRtnIt=ARIst+nRIt=A∧RIΔst+nRIt=A∧RIsRIt+nRIt,where(21)A^RI=T¯θ1,T¯θ2,…,T¯θK∈R2M×2FK,Δ=blkdiagΛθ1,φ1,G,Λθ2,φ2,G,…,ΛθK,φK,G,where A∧RI denotes the real-valued block extended array manifold matrix formed by T¯θkk=1,2,…,K. It separates angle parameter from disturbance factors, such as mutual coupling coefficients and noncircular phase, making it only depend on DOAs information. The block diagonal matrix Δ∈ℝ2FK×K is combined with the real signal vector st to construct a novel block signal vector sRIt. i.e., sRIt=Δst∈ℝ2FK×1. It can be found that the 2Fk−2F+1 th to 2Fk th rows in sRIt correspond to k th element in st. Besides, both the new extended array manifold matrix A∧RI and the corresponding signal vector sRIt in (20) have block structures for each target, unlike the original coupled extended steering matrix ARI and the signal vector st in (10).


## 3. Real-Valued Weighted Subspace Fitting with Block Sparse Recovery

### 3.1. The Subspace Fitting Framework with Optimal Weighted Matrix

It has been analyzed that the real-valued coupled extended array manifold matrixARI still spans the same range subspace as the signal subspace Es. Through combining the basic theory of linear algebra and the fact that Δ in (21) is a column full rank matrix, it can be deduced that the signal subspace Es is a subset of the range space of the real-valued block extended array manifold matrix A∧RI, which satisfies(22)Es=ARIU=A∧RIΔU=A∧RIU∧,where U∧=ΔU∈ℝ2FK×K is a block diagonal matrix with column full rank and consists of K subblocks U∧kk=1,2,…,K. As the k th subblock, U∧k is composed of the 2Fk−2F+1 th to 2Fk th rows in U∧ corresponding to k th row in U. Whereas, (22) will be invalid when there are disturbances such as noise in the array output.As revealed in [39], numerous prevalent works on DOA estimation can be viewed as the subspace fitting problem in a general sense. Then, a subspace fitting framework is given as follows:(23)θ∧,U↔=argminθ,U˜EsW1/2−A∧RIθU˜F2,where A∧RI is parameterized by θ. W∈ℝK×K represents a positive definite weighted matrix depending on the distinct calculation ways and affecting the asymptotic characteristics of fitting error. According to [39], it has been revealed that there exists an optimal weighted matrix that asymptotically minimizes the fitting error variance in the target directions and satisfies Wopt=Ωs−σ∧2IK2Ωs−1. σ∧2 denotes any consistent estimate of noise variance, which can be achieved by averaging 2M−K smaller eigenvalues of R. Highlighting that when W=Wopt, (23) describes the optimal subspace fitting issue, defined as weighted subspace fitting (WSF) problem [40].It is emphasized thatA∧RI and U˜ can be separated in the process of subspace fitting [38]. Meanwhile, the parameter we care about is A∧RI, not U˜. Therefore, the least square solution of U˜ can be solved by fixing A∧RI. i.e.,(24)U↔=A∧RI†θEsW1/2.Then bringing (24) back to (23), yields(25)θ∧=argminθtrPA∧RIθ⊥EsWoptEsH=argminθΥθ,where PA∧RIθ⊥=I2M−PA∧RIθ=I2M−A∧RIθA∧RI†θ.In order to deeply study the subspace fitting issue structured by (25) from the perspective of sparse reconstruction, the spatial domain is evenly discretized to form an over-complete dictionary A¯RI. A¯RI takes the following form:(26)A¯RI=T¯θ¯1,T¯θ¯2,…,T¯θ¯N∈ℝ2M×2FN,where θ¯=θ¯1,θ¯2,…,θ¯N represents a sampling grid point set and N indicates the number of grid points. It is noted that compared with M and K, N is sufficiently large in this paper, so that the grid-off issue is not considered here. Through combining (22) with Wopt, EsWopt1/2=A∧RIU∧Wopt1/2=A∧RIU˜ can be structured. However, it should be pointed out that such relationship is mathematically strict only under the condition of infinite snapshots. Then, based on the over-complete dictionary in (26), EsWopt1/2 can be sparsely represented as(27)EsWopt1/2=A¯RIU¯,where U¯=U¯θ¯1T,U¯θ¯2T,…,U¯θ¯NTT denotes a block sparse matrix, whose n th subblock U¯θ¯n is made up of the 2Fn−2F+1 th to 2Fn th rows of U¯. Furthermore, the subblocks corresponding to the desired DOAs in U¯ are equal to those in U˜, while the rest are zero. i.e.,(28)U¯θ¯n=U˜θk,θ¯n∈θ1,θ2,…,θK0,θ¯n∉θ1,θ2,…,θK,where n=1,2,…,N and k=1,2,…,K.According to (28), it is known that there are only K nonzero subblocks in U¯ due to the existence of K targets. Therefore, the DOA estimation issue can be transformed into a block sparse recovery problem, in which DOAs can be estimated by determining the positions of nonzero subblocks in U¯.It can be discovered that block sparse matrixU¯ is critical for direction finding, which can be reconstructed via minimizing l0-norm. Then, a l0-norm optimization scheme is formed as(29)minrl20s.t.EsWopt1/2=A¯RIU¯,where a column vector rl2=r1l2,r2l2,…,rNl2T is introduced to describe sparsity. rnl2 is the n th element in rl2 and corresponds to the n th subblock of U¯, which can be computed by the l2-norm of the 2Fn−2F+1 th to 2Fn th rows in U¯. That is to say, rnl2=∑a=2Fn−2F+12Fn∑b=1KU¯a,b2, in which U¯a,b represents the element located at the a th row and b th column of U¯. Evidently, the sparsity of vector rl2 is the same as that of the block sparse matrix U¯.In general,l0-norm penalty can accurately describe the solutions sparsity in the process of sparse recovery. Whereas, l0-norm penalty is a nondeterministic polynomial (NP)-hard and nonconvex combinatorial optimization problem, so it is mathematically intractable. In this way, l0-norm convex relaxes to l1-norm to solve the above problem. Moreover, considering the fitting error caused by finite snapshots, the l1-norm penalty framework is ultimately restructured as(30)minrl21s.t.EsWopt1/2−A¯RIU¯F≤ε,where the regularization parameter ε means the upper limit of the subspace fitting error, that is utilized to guarantee robust DOA estimation. Inspired by (25), it is known that the subspace fitting error is equal to Υθ. It has been derived that function 2T/σ∧2Υθ asymptotically follows chi-square distribution with 2K′2M−K degrees of freedom when θ refers to the true DOAs [40]. Thereby, Υθ≤ε with a high confidence interval 1−ρ¯ is calculated to determine parameter ε, in which ρ¯=0.001 is chosen in this paper.


### 3.2. Reweighted Block Sparse Recovery for DOA Estimation

Through the sparse recovery framework achieved byl1-norm constrained optimization in (30), DOA estimation can indeed be obtained. However, l1-norm penalty is only the convex approximation of l0-norm minimization, resulting in limited recovery accuracy. Specifically speaking, the penalty imposed on larger coefficients is heavier than that imposed on smaller coefficients in the l1-norm penalty framework, unlike the democratic l0-norm constraint. Then, a weighted matrix is introduced to enhance the solutions sparsity, where the weights can be determined by constructing the penalty factors. Therefore, l1-norm can approximate l0-norm as much as possible.Following the principle of MUSIC-like approach in [36, 38], the orthogonality between the real-valued block extended array manifold and the noise subspace can be utilized to formulate a novel real-valued NC MUSIC-like spectrum function. The spectrum function can be expressed as(31)ZMUSICθ=1detT¯HθEnEnHT¯θ.Inspired by the discretized sampling gird pointsθ¯=θ¯1,θ¯2,…,θ¯N, the orthogonality between the over-complete dictionary and its noise subspace can be exploited to structure the weights. First, the over-complete dictionary in (26) can be categorized into two groups: A¯RI=A¯RI1,A¯RI2. It is supposed that A¯RI1 is composed of K block steering matrices corresponding to the interested DOAs, while A¯RI2 is formed by residual N−K subblocks. Then based on (31), the initial weights can be represented as(32)d∧n=detT¯Hθ¯nEnEnHT¯θ¯nn=1,2,…,N,where the weight d∧n indicates the determinant value corresponding to θ¯n. Then, a weighted matrix can be structured as(33)D=diagd,where D∈ℝN×N denotes a weighted matrix that depends on the vector d. Furthermore, d relates to the initial weights d∧nn=1,2,…,N and takes the following form:(34)d=d1,d2=d¯1,d¯2,…,d¯N=d∧1,d∧2,…,d∧Nmaxd∧1,d∧2,…,d∧N.It can be concluded that if the number of snapshots is sufficiently large, the weights ind1 corresponding to the interested DOAs are more likely to be zero, which are smaller than those in d2. Through exploiting the weighted matrix, larger coefficients are preserved by smaller weights, while smaller coefficients close to zero are punished by larger weights. Therefore, no matter how large or small the coefficients are, they can be punished more democratically, behaving like the fair l0-norm penalty. Eventually, by embedding the weighted matrix D, the reweighted scheme based on l1-norm optimization can be constructed as(35)minDrl21,s.t.EsWopt1/2−A¯RIU¯F≤ε.Thanks to second order cone (SOC) programming packages, like CVX, the optimization problem given in (35) can be successfully addressed. In this way, DOAs can be effectively estimated by detecting the positions of nonzero values in the reconstructed sparse vector rl2.Up to now, an efficient real-valued weighted subspace fitting algorithm with block sparse recovery has been proposed for DOA estimation of strictly noncircular sources with unknown mutual coupling, which can be summarized as algorithm 1.Algorithm 1: Real-valued weighted subspace fitting algorithm with block sparse recovery.
(1)
Input: The actual received signal yt in (6);(2)
Extract the real and imaginary parts ofyt based on (7) and (8) to formulate the real-valued coupled extended array output yRIt in (10);(3)
Calculate the sampling covariance matrixR of yRIt by (13);(4)
Perform eigenvalue decomposition onR to acquire the signal subspace Es and the noise subspace En in (14);(5)
Construct the optimal weighted matrixWopt according to (14);(6)
Form the over-complete dictionaryA¯RI in (26) by sparsely representing T¯θk of A^RI in (21) to develop a sparse representation model in (27);(7)
Structure the real-valued NC MUSIC-like weighted matrixD adopting (33) to enhance the solutions sparsity;(8)
Design the reweighted regularized framework based on WSF principle in (35);(9)
Output: The reconstructed real-valued sparse vector rl2;(10)
Perform a 1-D spectrum search to find theK maximum values for DOA estimation.


## 3.1. The Subspace Fitting Framework with Optimal Weighted Matrix

It has been analyzed that the real-valued coupled extended array manifold matrixARI still spans the same range subspace as the signal subspace Es. Through combining the basic theory of linear algebra and the fact that Δ in (21) is a column full rank matrix, it can be deduced that the signal subspace Es is a subset of the range space of the real-valued block extended array manifold matrix A∧RI, which satisfies(22)Es=ARIU=A∧RIΔU=A∧RIU∧,where U∧=ΔU∈ℝ2FK×K is a block diagonal matrix with column full rank and consists of K subblocks U∧kk=1,2,…,K. As the k th subblock, U∧k is composed of the 2Fk−2F+1 th to 2Fk th rows in U∧ corresponding to k th row in U. Whereas, (22) will be invalid when there are disturbances such as noise in the array output.As revealed in [39], numerous prevalent works on DOA estimation can be viewed as the subspace fitting problem in a general sense. Then, a subspace fitting framework is given as follows:(23)θ∧,U↔=argminθ,U˜EsW1/2−A∧RIθU˜F2,where A∧RI is parameterized by θ. W∈ℝK×K represents a positive definite weighted matrix depending on the distinct calculation ways and affecting the asymptotic characteristics of fitting error. According to [39], it has been revealed that there exists an optimal weighted matrix that asymptotically minimizes the fitting error variance in the target directions and satisfies Wopt=Ωs−σ∧2IK2Ωs−1. σ∧2 denotes any consistent estimate of noise variance, which can be achieved by averaging 2M−K smaller eigenvalues of R. Highlighting that when W=Wopt, (23) describes the optimal subspace fitting issue, defined as weighted subspace fitting (WSF) problem [40].It is emphasized thatA∧RI and U˜ can be separated in the process of subspace fitting [38]. Meanwhile, the parameter we care about is A∧RI, not U˜. Therefore, the least square solution of U˜ can be solved by fixing A∧RI. i.e.,(24)U↔=A∧RI†θEsW1/2.Then bringing (24) back to (23), yields(25)θ∧=argminθtrPA∧RIθ⊥EsWoptEsH=argminθΥθ,where PA∧RIθ⊥=I2M−PA∧RIθ=I2M−A∧RIθA∧RI†θ.In order to deeply study the subspace fitting issue structured by (25) from the perspective of sparse reconstruction, the spatial domain is evenly discretized to form an over-complete dictionary A¯RI. A¯RI takes the following form:(26)A¯RI=T¯θ¯1,T¯θ¯2,…,T¯θ¯N∈ℝ2M×2FN,where θ¯=θ¯1,θ¯2,…,θ¯N represents a sampling grid point set and N indicates the number of grid points. It is noted that compared with M and K, N is sufficiently large in this paper, so that the grid-off issue is not considered here. Through combining (22) with Wopt, EsWopt1/2=A∧RIU∧Wopt1/2=A∧RIU˜ can be structured. However, it should be pointed out that such relationship is mathematically strict only under the condition of infinite snapshots. Then, based on the over-complete dictionary in (26), EsWopt1/2 can be sparsely represented as(27)EsWopt1/2=A¯RIU¯,where U¯=U¯θ¯1T,U¯θ¯2T,…,U¯θ¯NTT denotes a block sparse matrix, whose n th subblock U¯θ¯n is made up of the 2Fn−2F+1 th to 2Fn th rows of U¯. Furthermore, the subblocks corresponding to the desired DOAs in U¯ are equal to those in U˜, while the rest are zero. i.e.,(28)U¯θ¯n=U˜θk,θ¯n∈θ1,θ2,…,θK0,θ¯n∉θ1,θ2,…,θK,where n=1,2,…,N and k=1,2,…,K.According to (28), it is known that there are only K nonzero subblocks in U¯ due to the existence of K targets. Therefore, the DOA estimation issue can be transformed into a block sparse recovery problem, in which DOAs can be estimated by determining the positions of nonzero subblocks in U¯.It can be discovered that block sparse matrixU¯ is critical for direction finding, which can be reconstructed via minimizing l0-norm. Then, a l0-norm optimization scheme is formed as(29)minrl20s.t.EsWopt1/2=A¯RIU¯,where a column vector rl2=r1l2,r2l2,…,rNl2T is introduced to describe sparsity. rnl2 is the n th element in rl2 and corresponds to the n th subblock of U¯, which can be computed by the l2-norm of the 2Fn−2F+1 th to 2Fn th rows in U¯. That is to say, rnl2=∑a=2Fn−2F+12Fn∑b=1KU¯a,b2, in which U¯a,b represents the element located at the a th row and b th column of U¯. Evidently, the sparsity of vector rl2 is the same as that of the block sparse matrix U¯.In general,l0-norm penalty can accurately describe the solutions sparsity in the process of sparse recovery. Whereas, l0-norm penalty is a nondeterministic polynomial (NP)-hard and nonconvex combinatorial optimization problem, so it is mathematically intractable. In this way, l0-norm convex relaxes to l1-norm to solve the above problem. Moreover, considering the fitting error caused by finite snapshots, the l1-norm penalty framework is ultimately restructured as(30)minrl21s.t.EsWopt1/2−A¯RIU¯F≤ε,where the regularization parameter ε means the upper limit of the subspace fitting error, that is utilized to guarantee robust DOA estimation. Inspired by (25), it is known that the subspace fitting error is equal to Υθ. It has been derived that function 2T/σ∧2Υθ asymptotically follows chi-square distribution with 2K′2M−K degrees of freedom when θ refers to the true DOAs [40]. Thereby, Υθ≤ε with a high confidence interval 1−ρ¯ is calculated to determine parameter ε, in which ρ¯=0.001 is chosen in this paper.


## 3.2. Reweighted Block Sparse Recovery for DOA Estimation

Through the sparse recovery framework achieved byl1-norm constrained optimization in (30), DOA estimation can indeed be obtained. However, l1-norm penalty is only the convex approximation of l0-norm minimization, resulting in limited recovery accuracy. Specifically speaking, the penalty imposed on larger coefficients is heavier than that imposed on smaller coefficients in the l1-norm penalty framework, unlike the democratic l0-norm constraint. Then, a weighted matrix is introduced to enhance the solutions sparsity, where the weights can be determined by constructing the penalty factors. Therefore, l1-norm can approximate l0-norm as much as possible.Following the principle of MUSIC-like approach in [36, 38], the orthogonality between the real-valued block extended array manifold and the noise subspace can be utilized to formulate a novel real-valued NC MUSIC-like spectrum function. The spectrum function can be expressed as(31)ZMUSICθ=1detT¯HθEnEnHT¯θ.Inspired by the discretized sampling gird pointsθ¯=θ¯1,θ¯2,…,θ¯N, the orthogonality between the over-complete dictionary and its noise subspace can be exploited to structure the weights. First, the over-complete dictionary in (26) can be categorized into two groups: A¯RI=A¯RI1,A¯RI2. It is supposed that A¯RI1 is composed of K block steering matrices corresponding to the interested DOAs, while A¯RI2 is formed by residual N−K subblocks. Then based on (31), the initial weights can be represented as(32)d∧n=detT¯Hθ¯nEnEnHT¯θ¯nn=1,2,…,N,where the weight d∧n indicates the determinant value corresponding to θ¯n. Then, a weighted matrix can be structured as(33)D=diagd,where D∈ℝN×N denotes a weighted matrix that depends on the vector d. Furthermore, d relates to the initial weights d∧nn=1,2,…,N and takes the following form:(34)d=d1,d2=d¯1,d¯2,…,d¯N=d∧1,d∧2,…,d∧Nmaxd∧1,d∧2,…,d∧N.It can be concluded that if the number of snapshots is sufficiently large, the weights ind1 corresponding to the interested DOAs are more likely to be zero, which are smaller than those in d2. Through exploiting the weighted matrix, larger coefficients are preserved by smaller weights, while smaller coefficients close to zero are punished by larger weights. Therefore, no matter how large or small the coefficients are, they can be punished more democratically, behaving like the fair l0-norm penalty. Eventually, by embedding the weighted matrix D, the reweighted scheme based on l1-norm optimization can be constructed as(35)minDrl21,s.t.EsWopt1/2−A¯RIU¯F≤ε.Thanks to second order cone (SOC) programming packages, like CVX, the optimization problem given in (35) can be successfully addressed. In this way, DOAs can be effectively estimated by detecting the positions of nonzero values in the reconstructed sparse vector rl2.Up to now, an efficient real-valued weighted subspace fitting algorithm with block sparse recovery has been proposed for DOA estimation of strictly noncircular sources with unknown mutual coupling, which can be summarized as algorithm 1.Algorithm 1: Real-valued weighted subspace fitting algorithm with block sparse recovery.
(1)
Input: The actual received signal yt in (6);(2)
Extract the real and imaginary parts ofyt based on (7) and (8) to formulate the real-valued coupled extended array output yRIt in (10);(3)
Calculate the sampling covariance matrixR of yRIt by (13);(4)
Perform eigenvalue decomposition onR to acquire the signal subspace Es and the noise subspace En in (14);(5)
Construct the optimal weighted matrixWopt according to (14);(6)
Form the over-complete dictionaryA¯RI in (26) by sparsely representing T¯θk of A^RI in (21) to develop a sparse representation model in (27);(7)
Structure the real-valued NC MUSIC-like weighted matrixD adopting (33) to enhance the solutions sparsity;(8)
Design the reweighted regularized framework based on WSF principle in (35);(9)
Output: The reconstructed real-valued sparse vector rl2;(10)
Perform a 1-D spectrum search to find theK maximum values for DOA estimation.


## 4. Simulation and Analysis

In this section, plenty of simulations are implemented and the corresponding results are exhibited to demonstrate the superior performance of the proposed method.


### 4.1. Simulation Scene

To demonstrate the superiority of the proposed method, the reweighted BSR method in [33] (defined as ReBSR method) and the joint reweighted sparsity-inducing method based on WSF principle in [38] (defined as WSF method) are chosen as the comparison methods. Meanwhile inspired by the Cramer–Rao bound (CRB) for noncircular signals of MIMO radar in [22], the array CRB is redrived for noncircular sources under unknown mutual coupling in this paper. In addition, the root mean square error (RMSE) is used to evaluate the estimation performance of all algorithms, which can be achieved by(36)RMSE=1JK∑j=1J∑k=1Kθ↔j,k−θk2,where θk stands for the real DOAs of k th target and θ↔j,k is estimated by θk in the j th Monte Carlo experiment. K refers to the number of radiating sources. The number of Monte Carlo experiments is set as J=100.In what follows, it is assumed thatM=8 antennas form a ULA, each of them is separated by half-wavelength spacing. There are K=2 narrowband uncorrelated strictly noncircular sources incident on the ULA from different directions in the far field, whose DOAs are denoted as θ1=−10° and θ2=2°, respectively. Additionally, the mutual coupling matrix consists of H=3 nonzero coefficients including g1=0.6864−j0.0919 and g2=0.2079−j0.0603. The entire spatial domain from −90° to 90° is uniformly sampled at the grid interval of 0.1°.


### 4.2. Simulation Results

Figure2 gives the spatial spectra for all different methods, in which SNR is set to −5 dB and the number of snapshots is fixed at 100. According to Figure 2, it can be observed that all methods form spectral peaks at the true DOAs positions. That is to say, these algorithms can realize direction finding of noncircular sources in the case of unknown mutual coupling. Furthermore, ReBSR method has the least shape peaks, the highest side-lobe, and the farthest from the real DOAs, while the proposed approach has the sharpest peaks, the lowest side-lobe, and the closest to the desirable DOAs. In this way, the proposed method outperforms other approaches in terms of resolution and accuracy.Figure 2
The spatial spectra for all methods.Figures3 and 4 indicate RMSE versus SNR and PSD versus SNR for distinct algorithms, respectively. In Figure 3, the number of snapshots is chosen as T=100. On the one hand, as displayed in Figure 3, the RMSE of these methods gradually decreases as SNR increases. Therefore, these methods can achieve improved performance by enhancing the signal environment. Moreover, the main difference between the proposed method and WSF algorithm is whether to perform real-valued conversion, so its impact on the estimation accuracy may not be obvious. In other words, it is reasonable to assume that their performance is similar. At the same time, the RMSE of ReBSR algorithm is larger than that of other noncircular methods, which is mainly due to its inability to take advantage of noncircularity, unlike the other two noncircular algorithms. On the other hand, as given in Figure 4, PSD refers to the successful detection rate for all Monte Carlo running experiments. And when the error absolute value between the true DOA θk and the estimated DOA θ↔k is less than 0.3°, the target detection is considered successful. From Figure 4, within the selected SNR range, the PSD of the proposed method and WSF approach are much higher than that of ReBSR algorithm. Additionally, they can be the first to achieve 100% PSD compared to ReBSR method. In conclusion, the proposed method has advantage over ReBSR algorithm and similar performance to WSF approach.Figure 3
RMSE versus SNR for all methods.Figure 4
PSD versus SNR for all methods.Figure5 depicts RMSE versus snapshots for distinct methods, when SNR is fixed at SNR=0 dB. As shown in Figure 5, the overall simulation trend is similar to that of Figure 3. As the number of snapshots increases, the RMSE of all distinct approaches decreases. Furthermore, in terms of estimation accuracy, the proposed method is similar to WSF algorithm, better than ReBSR approach and closer to CRB.Figure 5
RMSE versus snapshots for all methods.Figure6 shows RMSE of the proposed method versus SNR for different number of antennas, in which T=100. From Figure 6, the RMSE is the largest when M=7, while the RMSE is the smallest when M=9, which means that if SNR is fixed, the RMSE of the proposed method decreases with the increase of the number of sensors. However, it is worth highlighting that the more sensors, the higher the estimation accuracy, the heavier the computational load, and even the stronger the antennas interaction. In other words, it is better to make a trade-off between effectiveness and efficiency.Figure 6
RMSE of the proposed method versus SNR for different number of antennas.Figure7 reveals the average simulation time required for the two noncircular methods versus gird interval, where SNR and snapshots are set to SNR=0 dB and T=100, respectively. In Figure 7, whether the proposed method or WSF algorithm, the larger the grid interval, the shorter the simulation time and the lower the computational burden. This is mainly because the number of sampling grid points decreases as the grid interval increases. In addition, it is evident that the proposed method requires much less time than WSF algorithm, which means that the proposed method is superior to WSF algorithm in simulation time, although they all belong to noncircular algorithms. The main reason is that the proposed method converts complex domain date into real domain data to speed up direction finding for meeting the real-time requirement as much as possible, different from WSF approach in the complex domain. In this way, it takes less time to efficiently achieve high-precision DOA estimation, which is more suitable for practical applications.Figure 7
Average simulation time versus grid interval.


## 4.1. Simulation Scene

To demonstrate the superiority of the proposed method, the reweighted BSR method in [33] (defined as ReBSR method) and the joint reweighted sparsity-inducing method based on WSF principle in [38] (defined as WSF method) are chosen as the comparison methods. Meanwhile inspired by the Cramer–Rao bound (CRB) for noncircular signals of MIMO radar in [22], the array CRB is redrived for noncircular sources under unknown mutual coupling in this paper. In addition, the root mean square error (RMSE) is used to evaluate the estimation performance of all algorithms, which can be achieved by(36)RMSE=1JK∑j=1J∑k=1Kθ↔j,k−θk2,where θk stands for the real DOAs of k th target and θ↔j,k is estimated by θk in the j th Monte Carlo experiment. K refers to the number of radiating sources. The number of Monte Carlo experiments is set as J=100.In what follows, it is assumed thatM=8 antennas form a ULA, each of them is separated by half-wavelength spacing. There are K=2 narrowband uncorrelated strictly noncircular sources incident on the ULA from different directions in the far field, whose DOAs are denoted as θ1=−10° and θ2=2°, respectively. Additionally, the mutual coupling matrix consists of H=3 nonzero coefficients including g1=0.6864−j0.0919 and g2=0.2079−j0.0603. The entire spatial domain from −90° to 90° is uniformly sampled at the grid interval of 0.1°.


## 4.2. Simulation Results

Figure2 gives the spatial spectra for all different methods, in which SNR is set to −5 dB and the number of snapshots is fixed at 100. According to Figure 2, it can be observed that all methods form spectral peaks at the true DOAs positions. That is to say, these algorithms can realize direction finding of noncircular sources in the case of unknown mutual coupling. Furthermore, ReBSR method has the least shape peaks, the highest side-lobe, and the farthest from the real DOAs, while the proposed approach has the sharpest peaks, the lowest side-lobe, and the closest to the desirable DOAs. In this way, the proposed method outperforms other approaches in terms of resolution and accuracy.Figure 2
The spatial spectra for all methods.Figures3 and 4 indicate RMSE versus SNR and PSD versus SNR for distinct algorithms, respectively. In Figure 3, the number of snapshots is chosen as T=100. On the one hand, as displayed in Figure 3, the RMSE of these methods gradually decreases as SNR increases. Therefore, these methods can achieve improved performance by enhancing the signal environment. Moreover, the main difference between the proposed method and WSF algorithm is whether to perform real-valued conversion, so its impact on the estimation accuracy may not be obvious. In other words, it is reasonable to assume that their performance is similar. At the same time, the RMSE of ReBSR algorithm is larger than that of other noncircular methods, which is mainly due to its inability to take advantage of noncircularity, unlike the other two noncircular algorithms. On the other hand, as given in Figure 4, PSD refers to the successful detection rate for all Monte Carlo running experiments. And when the error absolute value between the true DOA θk and the estimated DOA θ↔k is less than 0.3°, the target detection is considered successful. From Figure 4, within the selected SNR range, the PSD of the proposed method and WSF approach are much higher than that of ReBSR algorithm. Additionally, they can be the first to achieve 100% PSD compared to ReBSR method. In conclusion, the proposed method has advantage over ReBSR algorithm and similar performance to WSF approach.Figure 3
RMSE versus SNR for all methods.Figure 4
PSD versus SNR for all methods.Figure5 depicts RMSE versus snapshots for distinct methods, when SNR is fixed at SNR=0 dB. As shown in Figure 5, the overall simulation trend is similar to that of Figure 3. As the number of snapshots increases, the RMSE of all distinct approaches decreases. Furthermore, in terms of estimation accuracy, the proposed method is similar to WSF algorithm, better than ReBSR approach and closer to CRB.Figure 5
RMSE versus snapshots for all methods.Figure6 shows RMSE of the proposed method versus SNR for different number of antennas, in which T=100. From Figure 6, the RMSE is the largest when M=7, while the RMSE is the smallest when M=9, which means that if SNR is fixed, the RMSE of the proposed method decreases with the increase of the number of sensors. However, it is worth highlighting that the more sensors, the higher the estimation accuracy, the heavier the computational load, and even the stronger the antennas interaction. In other words, it is better to make a trade-off between effectiveness and efficiency.Figure 6
RMSE of the proposed method versus SNR for different number of antennas.Figure7 reveals the average simulation time required for the two noncircular methods versus gird interval, where SNR and snapshots are set to SNR=0 dB and T=100, respectively. In Figure 7, whether the proposed method or WSF algorithm, the larger the grid interval, the shorter the simulation time and the lower the computational burden. This is mainly because the number of sampling grid points decreases as the grid interval increases. In addition, it is evident that the proposed method requires much less time than WSF algorithm, which means that the proposed method is superior to WSF algorithm in simulation time, although they all belong to noncircular algorithms. The main reason is that the proposed method converts complex domain date into real domain data to speed up direction finding for meeting the real-time requirement as much as possible, different from WSF approach in the complex domain. In this way, it takes less time to efficiently achieve high-precision DOA estimation, which is more suitable for practical applications.Figure 7
Average simulation time versus grid interval.


## 5. Conclusion

In this paper, the scenario of strictly noncircular sources in the presence of unknown mutual coupling is concerned, and then a real-valued reweighted block sparse recovery framework achieved by WSF principle is structured for DOA estimation. In the proposed method, a real-valued coupled extended array output is first constructed by connecting the real and imaginary parts of the received data. Then, the real-valued block extended sparse recovery model is formed by exploiting the parameterized decoupling thought to avoid the influences of unknown mutual coupling and noncircular phase. Afterwards, a robust WSF approach is explored to recover the real-valued block sparse matrix for DOA estimation, where a real-valued NC MUSIC-like weighted matrix is further embedded to reinforce the solutions sparsity. Additionally, the upper bound of subspace fitting error is reported as well. Thanks to noncircularity, parameterized decoupling operation, and reweighted measure, the proposed method can not only provide desirable estimation accuracy, but also bear low computational burden. Extensive experiment results validate the effectiveness and efficiency of the proposed method for strictly noncircular sources with unknown mutual coupling.



---

*Source: 7199603-2023-01-12.xml*