# Blind Direction-of-Arrival Estimation with Uniform Circular Array in Presence of Mutual Coupling

Authors: Song Liu; Lisheng Yang; Shizhong Yang; Qingping Jiang; Haowei Wu
Journal: International Journal of Antennas and Propagation
 (2016)

---

## Abstract

A blind direction-of-arrival (DOA) estimation algorithm based on the estimation of signal parameters via rotational invariance techniques (ESPRIT) is proposed for a uniform circular array (UCA) when strong electromagnetic mutual coupling is present. First, an updated UCA model with mutual coupling in a discrete Fourier transform (DFT) beam space is deduced, and the new manifold matrix is equal to the product of a centrosymmetric diagonal matrix and a Vandermonde-structure matrix. Then we carry out blind DOA estimation through a modified ESPRIT method, thus avoiding the need for spatial angular searching. In addition, two mutual coupling parameter estimation methods are presented after the DOAs have been estimated. Simulation results show that the new algorithm is reliable and effective especially for closely spaced signals.


---

## Body

## 1. Introduction

Uniform circular array based DOA estimation methods are always attractive because the UCA has a special symmetric structure that provides almost the same resolution ability along the 360° azimuth angle domain. Except for conventional DOA estimation methods such as beam space searching or the Capon method [1], many other methods can be applied, such as the commonly used multiple signal classification (MUSIC) [2] or the well-known ESPRIT [3] in phase mode space. More effective methods such as UCA-RB-MUSIC and UCA-ESPRIT [4] have been introduced, and the mapping error reducing method was developed [5].However, the electromagnetic mutual coupling effect cannot be ignored in a real array. Generally, this effect will severely degrade the performance of the above methods [6]. A classic DOA estimation method based on an iterative search technique was proposed to estimate the DOA and the mutual coupling matrix (MCM) parameters jointly [7], but it has a high computation cost. The rank reduction (RARE) method introduced by [7] was further developed to obtain the DOA estimate for a UCA [8–10], but an angular search is still needed. Moreover, angular ambiguities are challenges for RARE-based blind methods. For example, two closely spaced signals cannot be differentiated according to the UCA-RARE [8] spectrum because there is a spurious peak at the middle of the two true angles.In this paper, we design a new modified ESPRIT method for UCA to estimate the azimuth angle when mutual coupling is present. The proposed method is blind to mutual coupling and can completely avoid angular searching. Reliable DOA estimates can be obtained, especially for closely spaced signals. In addition, we will introduce two methods to estimate the MCM parameters once the DOA values are calculated.


## 2. UCA Model with Mutual Coupling

SupposeM-element UCA has a radius r (Figure 1). All of the antenna elements are identical, and there are D far-field narrow signals impinging from θj,j=1,…,D which are the parameters to be estimated. The snapshot can be written as(1)zk=CAsk+nk,k=1,…,K.K is the total sampling number and A is the manifold matrix:(2)A=aθ1,…,aθD,aθjm=fθjexp⁡i2πrλcos⁡θj-2m-1πM,where m=1,…,M. λ is the wavelength of the signals, and aθj is the manifold vector for the jth signal. f(θj) is the common directivity factor of the antenna elements. Since all the antenna elements are identical, fθj can be normalized as 1, and we use effective SNR which has already included the directivity factor. sk is the signal sample vector:(3)sk=s1k,…,sDkT.C is the MCM, which is a symmetric circular Toeplitz matrix with the form(4)C=c0c1…cN…c1c1c0…cN-1…c2⋮⋮…⋮…⋮cNcN-1…c0…cN-1⋮⋮…⋮…⋮c1c2…cN-1…c0,N=M2.· is the flooring function and c0 is normalized as 1. For ideal dipole antenna array (Figure 1), induced EMF (Electromotive Force) method can give a close form of self-impedance and mutual impedance and thus MCM can be calculated according to Gupta and Ksienski’s formulation [6].Figure 1
Uniform circular array.Suppose the signals are uncorrelated andn(k) is white Gaussian noise:(5)Esk=0,EskskH=Rs,Enk=0,EnknkH=σw2I.


## 3. Algorithm

### 3.1. UCA Model in DFT Beam Space

First we introduce the DFT of the MCM for UCA.Lemma 1 (see [11–13]).
IfC is a circular matrix with its first column vector as c=[c0,c1,…,cM-1]T and F is a Fourier matrix(6)F=w0,w1,…,wM-1,(7)wj=1M1,exp⁡i2πjM,…,exp⁡i2πjM-1MT,where j=0,…,M-1, then Λ=FCFH is a diagonal matrix with entries as C’s eigenvalues:(8)μm=Λm+1m+1=∑j=0M-1cjexp⁡i2πjmM,in which m=0,1,…,M-1. For the Fourier matrix F(9)F=FT,FH=F∗.
IfC is also symmetric, this means that(10)C=CTorcj=cM-j,j=1,…,M-1.
Then we can rewriteΛ as(11)μm=wmTCwm∗=wmTCwm∗T=wmHCwm,0=wmTCwj∗m≠jand we have(12)μM-m=wM-mHCwM-m=wmTCwm∗=μm,where m=1,…,M-1. From (8) and (12), we get a linear equation which c′and μ′ should satisfy for even M:(13)Θc′=μ′,where(14)c′≔c0,…,cM/2T,μ′≔μ0,…,μM/2T,(15)Θ≔12⋯2112cos⁡2πM·1·1⋯2cos⁡2πM·1·M-22-112cos⁡2πM·2·1⋯2cos⁡2πM·2·M-221⋮⋮⋯⋮⋮12cos⁡2πM·M-22·1⋯2cos⁡2πM·M-22·M-22-112cos⁡2πM·M2·1⋯2cos⁡2πM·M2·M-221.
Therefore, if we determine the estimate ofμ′, then we can obtain the estimate of the mutual coupling parameters c′ by (13). In addition, there are similar equations where M is odd.
According to [4], we set(16)β=2πrλ,L=β,where L is the number of excited phase modes. We define another Fourier matrix F′, which is different from (6):(17)F′=F1,F2,(18)F1=w-L,…,wL,F2=wL+1,…,wM-L-1.
Then the snapshot in DFT beam space is(19)z~k=F′Hzk=F′HCF′F′HAsk+F′Hnk,if we write F′HCF′ as(20)F′HCF′=ΓC00ΓC′,where (see (11), (12), and (18))(21)ΓC=diag⁡μL,…,μ1,μ0,μ1,…,μL.
In addition, we can obtain [4](22)F′HA=F1HAF2HA=ΓJA~F2HAwith(23)ΓJ=Mdiag⁡i-LJ-Lβ,…,iLJLβ.Jl· is l order first-kind Bessel function. According to its property(24)J-lβ=-1lJlβwe get(25)ΓJ=Mdiag⁡αL,αL-1,…,α0,…,αL-1,αL,where(26)αl=ilJlβ,l=0,1,…,L.A~ is the updated manifold matrix with a Vandermonde structure(27)A~=exp⁡-iLθ1⋯exp⁡-iLθD⋮⋯⋮exp⁡iLθ1⋯exp⁡iLθD.
Finally, we obtain the snapshot in DFT beam space(28)z~k=ΓA~ΓC′F2HAsk+n~kin which(29)Γ=ΓCΓJ,En~kn~kH=σw2I.


### 3.2. DOA Estimation and MCM Calculation

We carry out the eigendecomposition on the covariance matrixRz~ and obtain the signal subspace U~s, which consists of eigenvectors corresponding to D maximum eigenvalues. Select the first L+1 rows of U~s as U~s1 and the second L+1 rows as U~s2. Select the first L+1 rows of A~ as A~1 and the second L+1 rows as A~2. Select the first L+1×(L+1) diagonal matrix of Γ as Γ1 and the second L+1×(L+1) diagonal matrix as Γ2. Use the same notations for other matrices ΓJ1, ΓJ2, ΓC1, and ΓC2. Then we have(30)U~s1=Γ1A~1T,U~s2=Γ2A~2T,where T is a nonsingular matrix. Thus we get(31)Γ12U~s2=U~s1Ψwith(32)Γ12=Γ1Γ2-1=ΓJ1ΓC1ΓJ2-1ΓC2-1=diag⁡d,d=μLαLμL-1αL-1,μL-1αL-1μL-2αL-2,…,μ1α1μ0α0,μ0α0μ1α1T,Ψ=T-1diagexp⁡iθ1,…,exp⁡iθDT.Since we use sampled data, we should replaceU~s with U~^s and define the object function as(33)d^,Ψ^=argmind,ΨdiagdU~^s2-U~^s1Ψs.t.d^Ld^L+1=1.We use the solution to the above equations from [14, 15]. Consider(34)d^=argminddHQds.t.d^Ld^L+1=1(35)Ψ^=U~^s1HU~^s1-1U~^s1Hdiag⁡d^U~^s2with(36)Q=PU~^s1⊥⊙U~^s2HU~^s2T,where PU~^s1⊥=I-U~^s1U~^s1HU~^s1-1U~^s1H and ⊙ is the Hadamard product. If ψmin is the eigenvector of Q corresponding to its minimum eigenvalue, then we have(37)d^=ψminρ,ρ2=ψminLψminL+1.Equation (37) indicates that there is a phase π ambiguity for ρ, and thus it will introduce a π ambiguity to the DOA values through the eigenvalues of Ψ^. However, this ambiguity can be cleared by comparing the RARE spectrum [9] values on θ^j and θ^j+π. We mark the estimated vector as d^opt without π ambiguity.Finally, we can determine the DOA estimates through the eigenvalues ofΨ^. Following is the detailed procedure of the blind method for DOA estimation of θj:(1)
Calculate the sample covariance matrixR^z and do eigendecomposition. Get the estimated signal subspace U^s and noise subspace U^n.(2)
Get the updated signal subspaceU~^s. Select the first L+1 rows as U~^s1 and the second L+1 rows as U~^s2:(38)U~^s=F′HU^s.(3)
Estimate the vectord^ and fitting matrix Ψ^ according to (34)–(36). Carry out eigendecomposition on Ψ^ and get the DOA estimates θ^j.(4)
Compare the blind RARE spectral values [9] on θ^j and θ^j+π and clear the ambiguity. Output the final DOA estimates.We label the above mentioned method as “Blind-m1-half” because we only selectL+1 rows of U~^s. In addition, we can select the first 2L and the second 2L rows from U~^s as U~^s1 and U~^s2. This method is labeled “Blind-m1-full.” Furthermore, we can select the first L and the third L rows from U~^s as U~^s1 and U~^s2 or select the first 2L-1 and the third 2L-1 rows from U~^s as U~^s1 and U~^s2. We label the two methods as “Blind-m2-half” and “Blind-m2-full,” respectively. We will carry out simulations on these four methods in Section 4.Now we can estimate the MCM parameters once we getθ^.Method 1.
Ifcj=0,j>L, then, according to (32) and (13), we get(39)diag⁡μ^Lμ^L-1,…,μ^0μ^1=ΓJ2ΓJ1-1diag⁡d^opt,(40)c^′=Θ^-1μ^′Θ^-1μ^′1,where c^′=[1,c^1,…,c^L]T, μ^′=[μ^0,μ^1,…,μ^L]T, and Θ^=Θ1:L+1,1:L+1), which means that Θ^ consists of the first L+1 rows and L+1 columns of Θ (see (15)). In (40), μ^0 is a nonzero parameter, but it will be cleared by the normalization operation.Method1 can only estimate L mutual coupling parameters.Method 2.
In addition, we can obtain a more accurate estimate ofc by the method introduced in [7]. Consider(41)c^′=argmincc′HΩc′s.t.c^′1=c0=1Ω=ΥHaθ^U^nU^nHΥaθ^,where c^′=[1,c^1,…,c^N]T. If ϕmin is the eigenvector of Ω corresponding to its minimum eigenvalue, then(42)c^′=ϕminϕmin1.
We give a simpler expression ofΥaθ^∈CM×(N+1) than that in [7, 9]. Consider(43)Υ:,1=aθ^,Υ:,j=circshiftaθ^,j-1+circshiftaθ^,1-j,j=2,3,…,N,Υ:,N+1=circshiftaθ^,N,for  evenNcircshiftaθ^,N+circshiftaθ^,-N,for  oddN,where circshiftaθ^,j is a function that shifts aθ^ by j times circularly. If j<0, then it shifts aθ^ in the counter direction.


## 3.1. UCA Model in DFT Beam Space

First we introduce the DFT of the MCM for UCA.Lemma 1 (see [11–13]).
IfC is a circular matrix with its first column vector as c=[c0,c1,…,cM-1]T and F is a Fourier matrix(6)F=w0,w1,…,wM-1,(7)wj=1M1,exp⁡i2πjM,…,exp⁡i2πjM-1MT,where j=0,…,M-1, then Λ=FCFH is a diagonal matrix with entries as C’s eigenvalues:(8)μm=Λm+1m+1=∑j=0M-1cjexp⁡i2πjmM,in which m=0,1,…,M-1. For the Fourier matrix F(9)F=FT,FH=F∗.
IfC is also symmetric, this means that(10)C=CTorcj=cM-j,j=1,…,M-1.
Then we can rewriteΛ as(11)μm=wmTCwm∗=wmTCwm∗T=wmHCwm,0=wmTCwj∗m≠jand we have(12)μM-m=wM-mHCwM-m=wmTCwm∗=μm,where m=1,…,M-1. From (8) and (12), we get a linear equation which c′and μ′ should satisfy for even M:(13)Θc′=μ′,where(14)c′≔c0,…,cM/2T,μ′≔μ0,…,μM/2T,(15)Θ≔12⋯2112cos⁡2πM·1·1⋯2cos⁡2πM·1·M-22-112cos⁡2πM·2·1⋯2cos⁡2πM·2·M-221⋮⋮⋯⋮⋮12cos⁡2πM·M-22·1⋯2cos⁡2πM·M-22·M-22-112cos⁡2πM·M2·1⋯2cos⁡2πM·M2·M-221.
Therefore, if we determine the estimate ofμ′, then we can obtain the estimate of the mutual coupling parameters c′ by (13). In addition, there are similar equations where M is odd.
According to [4], we set(16)β=2πrλ,L=β,where L is the number of excited phase modes. We define another Fourier matrix F′, which is different from (6):(17)F′=F1,F2,(18)F1=w-L,…,wL,F2=wL+1,…,wM-L-1.
Then the snapshot in DFT beam space is(19)z~k=F′Hzk=F′HCF′F′HAsk+F′Hnk,if we write F′HCF′ as(20)F′HCF′=ΓC00ΓC′,where (see (11), (12), and (18))(21)ΓC=diag⁡μL,…,μ1,μ0,μ1,…,μL.
In addition, we can obtain [4](22)F′HA=F1HAF2HA=ΓJA~F2HAwith(23)ΓJ=Mdiag⁡i-LJ-Lβ,…,iLJLβ.Jl· is l order first-kind Bessel function. According to its property(24)J-lβ=-1lJlβwe get(25)ΓJ=Mdiag⁡αL,αL-1,…,α0,…,αL-1,αL,where(26)αl=ilJlβ,l=0,1,…,L.A~ is the updated manifold matrix with a Vandermonde structure(27)A~=exp⁡-iLθ1⋯exp⁡-iLθD⋮⋯⋮exp⁡iLθ1⋯exp⁡iLθD.
Finally, we obtain the snapshot in DFT beam space(28)z~k=ΓA~ΓC′F2HAsk+n~kin which(29)Γ=ΓCΓJ,En~kn~kH=σw2I.


## 3.2. DOA Estimation and MCM Calculation

We carry out the eigendecomposition on the covariance matrixRz~ and obtain the signal subspace U~s, which consists of eigenvectors corresponding to D maximum eigenvalues. Select the first L+1 rows of U~s as U~s1 and the second L+1 rows as U~s2. Select the first L+1 rows of A~ as A~1 and the second L+1 rows as A~2. Select the first L+1×(L+1) diagonal matrix of Γ as Γ1 and the second L+1×(L+1) diagonal matrix as Γ2. Use the same notations for other matrices ΓJ1, ΓJ2, ΓC1, and ΓC2. Then we have(30)U~s1=Γ1A~1T,U~s2=Γ2A~2T,where T is a nonsingular matrix. Thus we get(31)Γ12U~s2=U~s1Ψwith(32)Γ12=Γ1Γ2-1=ΓJ1ΓC1ΓJ2-1ΓC2-1=diag⁡d,d=μLαLμL-1αL-1,μL-1αL-1μL-2αL-2,…,μ1α1μ0α0,μ0α0μ1α1T,Ψ=T-1diagexp⁡iθ1,…,exp⁡iθDT.Since we use sampled data, we should replaceU~s with U~^s and define the object function as(33)d^,Ψ^=argmind,ΨdiagdU~^s2-U~^s1Ψs.t.d^Ld^L+1=1.We use the solution to the above equations from [14, 15]. Consider(34)d^=argminddHQds.t.d^Ld^L+1=1(35)Ψ^=U~^s1HU~^s1-1U~^s1Hdiag⁡d^U~^s2with(36)Q=PU~^s1⊥⊙U~^s2HU~^s2T,where PU~^s1⊥=I-U~^s1U~^s1HU~^s1-1U~^s1H and ⊙ is the Hadamard product. If ψmin is the eigenvector of Q corresponding to its minimum eigenvalue, then we have(37)d^=ψminρ,ρ2=ψminLψminL+1.Equation (37) indicates that there is a phase π ambiguity for ρ, and thus it will introduce a π ambiguity to the DOA values through the eigenvalues of Ψ^. However, this ambiguity can be cleared by comparing the RARE spectrum [9] values on θ^j and θ^j+π. We mark the estimated vector as d^opt without π ambiguity.Finally, we can determine the DOA estimates through the eigenvalues ofΨ^. Following is the detailed procedure of the blind method for DOA estimation of θj:(1)
Calculate the sample covariance matrixR^z and do eigendecomposition. Get the estimated signal subspace U^s and noise subspace U^n.(2)
Get the updated signal subspaceU~^s. Select the first L+1 rows as U~^s1 and the second L+1 rows as U~^s2:(38)U~^s=F′HU^s.(3)
Estimate the vectord^ and fitting matrix Ψ^ according to (34)–(36). Carry out eigendecomposition on Ψ^ and get the DOA estimates θ^j.(4)
Compare the blind RARE spectral values [9] on θ^j and θ^j+π and clear the ambiguity. Output the final DOA estimates.We label the above mentioned method as “Blind-m1-half” because we only selectL+1 rows of U~^s. In addition, we can select the first 2L and the second 2L rows from U~^s as U~^s1 and U~^s2. This method is labeled “Blind-m1-full.” Furthermore, we can select the first L and the third L rows from U~^s as U~^s1 and U~^s2 or select the first 2L-1 and the third 2L-1 rows from U~^s as U~^s1 and U~^s2. We label the two methods as “Blind-m2-half” and “Blind-m2-full,” respectively. We will carry out simulations on these four methods in Section 4.Now we can estimate the MCM parameters once we getθ^.Method 1.
Ifcj=0,j>L, then, according to (32) and (13), we get(39)diag⁡μ^Lμ^L-1,…,μ^0μ^1=ΓJ2ΓJ1-1diag⁡d^opt,(40)c^′=Θ^-1μ^′Θ^-1μ^′1,where c^′=[1,c^1,…,c^L]T, μ^′=[μ^0,μ^1,…,μ^L]T, and Θ^=Θ1:L+1,1:L+1), which means that Θ^ consists of the first L+1 rows and L+1 columns of Θ (see (15)). In (40), μ^0 is a nonzero parameter, but it will be cleared by the normalization operation.Method1 can only estimate L mutual coupling parameters.Method 2.
In addition, we can obtain a more accurate estimate ofc by the method introduced in [7]. Consider(41)c^′=argmincc′HΩc′s.t.c^′1=c0=1Ω=ΥHaθ^U^nU^nHΥaθ^,where c^′=[1,c^1,…,c^N]T. If ϕmin is the eigenvector of Ω corresponding to its minimum eigenvalue, then(42)c^′=ϕminϕmin1.
We give a simpler expression ofΥaθ^∈CM×(N+1) than that in [7, 9]. Consider(43)Υ:,1=aθ^,Υ:,j=circshiftaθ^,j-1+circshiftaθ^,1-j,j=2,3,…,N,Υ:,N+1=circshiftaθ^,N,for  evenNcircshiftaθ^,N+circshiftaθ^,-N,for  oddN,where circshiftaθ^,j is a function that shifts aθ^ by j times circularly. If j<0, then it shifts aθ^ in the counter direction.


## 4. Simulations

Consider a 16-element half-wave dipole antenna UCA withf0=1.032GHz and r=0.7λ. Set L=4. The mutual coupling vector c is listed in Table 1. The theoretical parameters are calculated according to Gupta and Ksienski’s formulation [6], and small values are treated as cj=0,j=3,…,8. Real parameters should be estimated by array calibration method.Table 1
Mutual coupling parameters.
c
0
c
1
c
2
c
j, j=3,…,8
1
-
0.3204
+
0.1878
i
0.1496
+
0.1153
i
0Two signals are impinging from 35° and 45° with the same signal noise ratio (SNR). We apply the classic MUSIC, blind RARE [9], blind R-RARE [10], and UCA-RARE [8] methods to obtain azimuth estimates. The results are shown in Figure 2.Figure 2
MUSIC and RARE spectrum (SNR = 20 dB,K=1000, one trial).It shows that the MUSIC method and RARE-based blind methods cannot differentiate these two signals because RARE-based methods will introduce spurious estimates. We should obtain the initial estimates from the MUSIC spectrum to start the iterative method [7], but it is difficult to find two different spectral peaks from the spectrum of “MUSIC without MCM.”We apply the proposed blind method to the above example. The estimated average DOA absolute bias and root mean square error (RMSE) versus SNR are illustrated in Figures3–6 (K=10000, 100 trials). Method in [7] (ϵ=0.01) and the Cramér-Rao bound (CRB) with a known MCM are also presented [16].Figure 3
Average DOA absolute bias versus SNR for signal 1.Figure 4
Average DOA absolute bias versus SNR for signal 2.Figure 5
RMSE versus SNR for signal 1.Figure 6
RMSE versus SNR for signal 2.This shows that all the four proposed methods can give satisfactory estimates and that the tendency of the RMSE is the same as the CRB with an increase of SNR. Method Blind-m1-full and method Blind-m2-full are more effective than the other two methods because more rows ofU~^s are involved when we calculate the fitting matrix Ψ^ (see (35)). A comparison of simulations versus sampling number will give similar results. Method in [7] gives a biased estimate.Table2 lists the MCM parameter estimates based on the two methods introduced in Section 3.Table 2
Mutual coupling estimates (SNR = 30 dB,K=1000).
c
1
c
2
True value
-
0.3204
+
0.1878
i
0.1496
+
0.1153
i
Method1
-
0.3114
+
0.2053
i
0.1451
+
0.1179
i
Method2
-
0.3203
+
0.1877
i
0.1496
+
0.1154
iThis shows that Method2 can give more accurate results.


## 5. Conclusion

In DFT beam space, we utilize a modified ESPRIT algorithm to obtain a reliable DOA estimate when there is a severe mutual coupling effect. The new blind method is efficient because it avoids searching for the spectral peaks. For closely spaced signals, neither the classic MUSIC nor RARE-based methods provide a good estimate, while the proposed new method can produce an accurate estimate. Moreover, these direction estimates can be used for further MCM parameter estimation.



---

*Source: 8109013-2016-02-07.xml*