Unitary transformations for spherical harmonics MUSIC

Signal Processing 131 (2017) 441–446

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Unitary transformations for spherical harmonics MUSIC Qinghua Huang n, Guangfei Zhang, Longfei Xiang, Yong Fang Key Laboratory of Specialty Fiber Optics and Optical Access Networks, Shanghai University, Shanghai 200072, China

art ic l e i nf o

a b s t r a c t

Article history: Received 17 May 2016 Received in revised form 3 August 2016 Accepted 2 September 2016 Available online 8 September 2016

Spherical arrays have been widely used in direction-of-arrival (DOA) estimation in recent years. In this paper, two unitary transformations for DOA estimation using spherical arrays are developed to transform the complex search vector to a real vector. Furthermore, the covariance matrix is transformed into a real matrix which has real eigenvalues and eigenvectors. Therefore, we can develop a real spherical harmonics multiple signal classification (RSHMUSIC) method to estimate the locations of signals. By exploiting the phase shift characteristic of the real search vector, the computation cost can be further reduced. Compared with the complex spherical harmonics MUSIC, the proposed methods have a lower amount of computations and better performance. Simulation results demonstrate the performance of the real SHMUSIC method. & 2016 Elsevier B.V. All rights reserved.

Keywords: Spherical array Spherical harmonics Unitary transformation Conjugate symmetric property Multiple signal classification (MUSIC)

1. Introduction Multiple signal classification (MUSIC) [1] is a popular subspacebased technique for estimating direction-of-arrival (DOA) of incident signals. It decomposes the covariance matrix into signal subspace and noise subspace. Based on the noise subspace and the search vector, the MUSIC spectrum is formed to search the direction space. However, in array signal processing, the search vector is complex and the covariance matrix is a complex Hermitian matrix. Therefore, the eigendecomposition and the calculation of MUSIC search spectrum require complex computations. They need a large amount of computations. In order to reduce the complexity, a unitary transformation method [2,3] was proposed to transform the complex covariance matrix and the complex search vector into a real symmetric matrix and a real vector, respectively. It saves a significant amount of computations. However, these methods only suit an equally spaced linear array. The unitary transformation method was extended for a uniform rectangular array (URA) [4], two parallel uniform linear arrays (TP-ULAs) [5], and a uniform circular array (UCA) [6], respectively. A spherical array has been widely used in DOA estimation [7–10]. It has a symmetrical structure and can acquire high azimuth and elevation estimation accuracy. Complexvalued MUSIC method has been used to DOA estimation using spherical array [7,8]. It has a large amount of calculations. Kumar et al. proposed the spherical harmonics root-MUSIC [10]. They rewrote the steering vector to illustrate the Vandermonde n

Corresponding author. E-mail address: [email protected] (Q. Huang).

http://dx.doi.org/10.1016/j.sigpro.2016.09.002 0165-1684/& 2016 Elsevier B.V. All rights reserved.

structure of the array manifold in the spherical harmonics domain. However, this method is only able to get azimuth estimation with the determined elevation. Yan et al. proposed a real-valued MUSIC method [11] for arbitrary linear array geometries because the steering matrix satisfies A*(θ) = A( − θ), where θ is the signal directions. But the search vector of spherical arrays in the spherical harmonics domain does not have the conjugate property like that of linear arrays, so this method cannot be used for DOA estimation in the spherical harmonics domain. Each element of the real-valued beamspace manifold in [12] contains a sum term in which the degree changes from 0 to the highest order. Therefore, it does not have the independent phase shift property for azimuth or elevation. The computation cost cannot be further reduced. Yan et al. proposed the real-valued time-domain implementation of the broadband modal beamformer in the spherical harmonics domain [13]. The modal transformation is real-valued because the product of the spherical harmonic and its conjugate is real-valued according to the character of spherical harmonics. However, there is no spherical harmonics product in the array signal model (in the following model (6)) for DOA estimation. Therefore, utilizing the property of the steering vector, we want to construct unitary matrices to directly transform the complex steering vector of the model into a real one and then get a real-valued signal model. In this paper, we develop two unitary transformations for MUSIC using spherical arrays. The array model is firstly constructed in the spherical harmonics domain. The search vector depends on spherical harmonics. Therefore, the conjugate of the element in the search vector has a similar symmetric property as the corresponding spherical harmonic. This can be adopted to find a unitary transformation to make the search vector have a block conjugate symmetric property. Based on this property, the second

442

Q. Huang et al. / Signal Processing 131 (2017) 441–446

transformation are developed to transform the complex and conjugate symmetric search vector into a real vector. Furthermore, the covariance matrix can be transformed into a real matrix which has real eigenvalues and eigenvectors. Therefore, a real spherical harmonics MUSIC (RSHMUSIC) estimator can be developed to realize the spectral search only using real computations. This decreases the computational cost considerably.

2. Array model in spherical harmonics domain Let us consider an L-sensor spherical array with the radius R impinged by D narrowband far-field uncorrelated source signals in free field environment. The sensors are omnidirectional and no mutual coupling between them. The elevation and azimuth of the lth sensor are represented by ϕl , respectively, θl and Ωl = (θl, ϕl ) , l = 1, … , L . The dth source signal sd(t ) with the wavenumber k impinges on the array from the direction Φd = (ϑd, φd ), where ϑd and φd are the incident elevation and azimuth, respectively, d = 1, … , D. The array model can be described as D

x(t ) =

∑ a(Φd)sd(t ) + v(t ),

(1)

d=1 T

where is the array output vector, x(t ) = [x1(t ) , … , xL(t )] v(t ) = [v1(t ) , … , vL(t )]T is the additive white Gaussian noise vector, t is the snapshot index, and a(Φd ) is the dth steering vector which can be represented in spherical harmonics [12,14] as

a(Φd ) = Y(Ω)B(k )y(Φd ),

(2)

P (Φ) =

(2n + 1) (n − m)! m Pn (cos θ )eimϕ , 4π (n + m)!

(3)

(4)

(5)

Using the N-order inverse spherical Fourier transform (SFT) approximation, the array output vector x(t ) can be denoted as x(t ) = Y(Ω )x nm(t ), where x nm(t ) is the array output SFT coefficient vector [15,16]. In the same way, the noise signal can be expressed as v(t ) = Y(Ω )vnm(t ), where vnm(t ) is the noise SFT coefficient vector. By combining (1), (2), the inverse SFT representations for x(t ) and v(t ), and using the least squares criterion, we get the model in the spherical harmonics domain as follows D

∑ anm(Φd)sd(t ) + vnm(t ), d=1

(7)

3. Unitary transformations for SHMUSIC 3.1. Block centrohermitian matrix According to [17], a b × c complex centrohermitian matrix W satisfies W = Jb W*Jc , where Jb represents the b × bexchange matrix with ones on its anti-diagonal positions and zeroes elsewhere. An M × M matrix C can be divided into Γ blocks according to the row index. Each block has pτ rows, where τ = 1, 2, … , Γ and Γ

∑τ = 1 pτ = M . With the same division rule, the matrix can be also divided into Γ blocks according to the column index. Assuming each block of C is a centrohermitian matrix, a block diagonal exchange matrix can be defined as

¯J = blkdiag (J , J , …, J ). p p p 2

(8)

Γ

C = ¯JC*¯J .

(6)

where a nm(Φd ) = [a00(Φd ) , a1−1(Φd ) , a10(Φd ) , a11(Φd ) , … , aNN (Φd )]T is the new steering vector in the spherical harmonics domain with anm(Φd ) = bn(kR )[Ynm(Φd )]*.

(9)

For odd pτ , define a submatrix as

Ξ pτ

where θ ) is the associated Legendre polynomial, B(k ) = diag{b0(kR ) , b1(kR ) , b1(kR ) , b1(kR ) , … , bN (kR )}, bn(kR ) = 4πi njn (kR ) for open spherical arrays, i = −1 , jn (kR ) is the spherical Bessel function with the order n, and y(Φd ) is a spherical harmonic vector as follows:

x nm(t ) =

,

where blkdiag(∙) denotes the block diagonal matrix. The block diagonal exchange matrix ¯J satisfies ¯¯ JJ = IM , where IM is an M × M identity matrix. Definition 1: An M × M matrix C is called block centrohermitian if

Pnm(cos

y(Φd ) = [Y 00(Φd ), Y1−1(Φd ) , Y10(Φd ), Y11(Φd ) , ⋯ , YNN (Φd )]H .

2

e uH a nm(Φ)

where eu is the eigenvector corresponding to the uth eigenvalue H ( λ1 ≥ ⋯ ≥ λU ) of the covariance matrix R = E[x nm(t )x nm (t )] and a nm(Φ ) is the search vector [7]. Calculating the MUSIC spectrum requires complex computations because both the eigenvectors and the search vector are complex. In this paper, we develop two unitary transformations to change complex computations into real computations.

where N is the highest order of the spherical harmonics, U = (N + 1)2, Ynm(Ω ) is the spherical harmonic with order n and degree m defined as [15]

Y nm(Ω) = Y nm(θ , ϕ) =

1 U ∑u = D + 1

1

where Y(Ω ) is an L × U spherical harmonics matrix as follows

⎡ Y 0(Ω )Y −1(Ω )Y 0(Ω )Y1(Ω )⋯Y N (Ω ) ⎤ N 1 1 1 1 1 1 ⎥ ⎢ 0 1 1 ⎢ Y 0(Ω )Y −1(Ω )Y 0(Ω )Y1(Ω )⋯Y N (Ω )⎥ N 2 1 2 1 2 2 , Y(Ω) = ⎢ 0 2 1 ⎥ ⎥ ⎢⋮⋮⋮⋮⋱⋮ ⎢ Y 0(Ω )Y −1(Ω )Y 0(Ω )Y1(Ω )⋯Y N (Ω ) ⎥ ⎣ 0 L 1 L 1 L 1 L N L ⎦

The spherical harmonics MUSIC (SHMUSIC) spectrum can be expressed as

⎡ I 0(pτ − 1) /2 J(p − 1) /2 ⎤ τ ⎥ ⎢ (pτ − 1) /2 1 ⎢ T = 0(pτ − 1) /2 2 0(Tpτ − 1) /2 ⎥, ⎥ 2⎢ ⎢ −i J 0(pτ − 1) /2 iI(pτ − 1) /2 ⎥⎦ (pτ − 1) /2 ⎣

(10)

where 0(p − 1) /2 is a column vector containing (pτ − 1) /2zeros, and τ for even pτ , we have

Ξ pτ =

⎡ ⎤ 1 ⎢ I pτ /2 J pτ /2 ⎥ . 2 ⎢⎣ −i J p /2 iI pτ /2 ⎥⎦

(11)

τ

Considering all Γ blocks, we can construct a block diagonal matrix as

Ξ = blkdiag (Ξ p1, Ξ p2 , …, Ξ pΓ ).

(12) −1

It is easily shown that Ξ is unitary, i.e., Ξ

H

= Ξ , and satisfies

Ξ*¯J = Ξ.

(13)

Using the unitary matrix Ξ , we have the following theorem. Table 1 The block division of the estimated covariance matrix according to the row index. Order n Block index Row index

0 1 1

1 2 2–4

2 3 5–9

3 4 10–16

4 5 17–25

Q. Huang et al. / Signal Processing 131 (2017) 441–446

Theorem: For any M × M block centrohermitian matrix ΞCΞH is real. The proof of the theorem is in the Appendix.

C,

443

Table 2 The relationship between n and m of observation SFT coefficients and search vectors. in three LC-RSHMUSIC methods. LC-RSHMUSIC

g ¼13

g ¼ 15

g¼9

n¼ 0 n¼ 1 n¼ 2 n¼ 3 n¼ 4

m¼ 0 m¼ 0 m¼ 0, 7 2 m¼ 0, 7 2 m¼ 0, 7 2, 7 4

m¼ 0 m¼ 71 m¼ 0, 72 m¼ 71, 7 3 m¼ 0, 7 2, 7 4

m¼ 0

3.2. Two unitary transformations for SHMUSIC The search vector in (7) depends on spherical harmonics. According to the definition in (4), the spherical harmonic has the property of [Ynm(Φ )]* = ( − 1)mYn−m(Φ ) [15]. The element of the search vector has the similar property as

[a nm(Φ)]* = ( − 1)m + na n−m(Φ).

m¼ 0, 7 2 m¼ 0, 72, 7 4

(14)

Firstly we construct a unitary matrix to make the search vector have the block conjugate symmetric property. Each block of the unitary matrix can be constructed as follows:

⎡ In ⎤ ⎢ ⎥ n Qn = ⎢ i ⎥, ⎢⎣ Gn ⎥⎦

(15)

where Gn is an n × n diagonal matrix, G0 = 1, when n ≥ 1, Gn = diag{( − 1)n + 1, ( − 1)n + 2 , … , ( − 1)2n} . Then we can construct a block diagonal matrix as

Q = blkdiag (Q 0, Q 1, …, Q N).

(16) −1

H

Q is a unitary matrix, i.e., Q = Q . Multiplying (6) using the first unitary matrix Q from the left, we get D

x¯ nm(t ) =

∑ a¯ nm(Φd)sd(t ) + v¯nm(t ),

(17)

d=1

where x¯ nm(t ) = Qx nm(t ), a¯ nm(Φ ) = Qa nm(Φ ), and v¯nm(t ) = Qvnm(t ). For an order n, the corresponding subvector of the search vector is a n(Φ ) = [an−n(Φ ) , …, an−1(Φ ) , an0(Φ ) , an1(Φ ) , …, ann(Φ )]T . According to the definition in (15) of the nth block unitary matrix Q n , .

Q na n(Φ )

= [an−n(Φ ) , … , an−1(Φ ) , i nan0(Φ ) , ( − 1)n + 1an1(Φ ) , …, ( − 1)2nann(Φ )]T Therefore, the element a¯ nm(Φ ) of the new search vector a¯ nm(Φ ) can be expressed as

a¯ nm(Φ)

⎧ 4πi njn (kR)[Y nm(Φ)]* , m < 0, ⎪ ⎪ m = 0, ( − 1)n4πjn (kR)Y nm(Φ), =⎨ ⎪ n+ m n m ⎪ 4πi jn (kR)[Y n (Φ)]* , m > 0, ⎩ ( − 1) [a¯ nm(Φ )]*

Fig. 1. MUSIC spectrum of two signals with the same elevation versus azimuth.

1 ^ ^* R˜ = (R nm + JR nmJ). 2

(21)

Each block of R˜ is a centrohermitian submatrix [18]. Considering all the blocks, we have R˜ = JR˜ *J. So R˜ is a block centrohermitian matrix. As the dimension of each block in R˜ is odd, according to (10), we can construct a block diagonal matrix as

(18) a¯ n−m(Φ ).

= with the conjugate symmetric property For an order n, the degree m changes from –n to n. Therefore, the new search vector has the block conjugate symmetric property as

[a¯ nm(Φ)]* = Ja¯ nm(Φ),

(19)

where J is a block diagonal exchange matrix defined as

J = blkdiag (J1 , J3 , …, J2N + 1).

(20)

The covariance matrix is commonly obtained using a forward T ^ H ¯ nm(t )x¯ nm only averaged estimate R (t ) /T which is only nm = ∑t = 1 x Hermitian. Because the order of spherical harmonics is from 0 to ^ N, R nm dependent on the order can be divided into N þ1 blocks according to the row index. For n = 0, … , N , the (nþ 1)th block consists of the elements form the (n2 þ1)th rows to the (n þ1)2 th ^ rows. Take an example of N = 4 and divide R nm into 5 blocks ac^ cording to the row index shown in Table 1. In the same way, R nm can be divided into N þ 1 blocks according to the column index. Therefore, we use the block diagonal exchange matrix J to realize the block forward/backward (FB) averaged estimate of the covariance matrix

Fig. 2. MUSIC spectrum of two signals with the same azimuth versus elevation.

444

Q. Huang et al. / Signal Processing 131 (2017) 441–446

Fig. 4. Probability of correct elevation estimation versus SNR.

Fig. 3. Probability of correct azimuth estimation versus SNR.

Q¯ = blkdiag (Ξ1, Ξ3, …, Ξ2N + 1).

˜ ¯ H , the relations between the eigenvalues and eigenvectors of Q¯ RQ eigenvectors are λ¯u = λ¯′u and e′u = Q¯ e¯ u . The MUSIC spectrum in (7) becomes

P (Φ) =

=

U ∑u = D + 1 e¯ uH a¯ nm(Φ)

2

1 U ∑u = D + 1

e¯′uT a¯ ′nm(Φ)

2

=

1 U ∑u = D + 1 e¯′uT Q¯ a¯ nm(Φ)

2

,

RM SE of az im uth (degree)

Because the FB averaged estimate R˜ is a block centrohermitian matrix, it can be transformed into a real matrix by the second unitary transformation matrix Q¯ using the above theorem. Fur˜ ¯ H is symmetric since R˜ is Hermitian. Thus, Q¯ RQ ˜ ¯ H is thermore, Q¯ RQ real and symmetric. Let λ¯u and e¯ u ( u = 1, … , U ) be the eigenvalues and eigenvectors of R˜ , λ¯′u and e′u be the real eigenvalues and real

1

5

(22)

3

2

1

0

0

5

10

15

20

SNR (dB) Fig. 5. RMSE of azimuth estimation versus SNR.

(24)

where Re{} and Im{} denotes the real part and imaginary part of a complex value, respectively. From (24), we can see that the element of our proposed real search vector contains the real part or imaginary part of the spherical harmonics dependent on the single degree m and avoids the summation of all terms related to positive m in the real beamspace manifold of [12]. In computing the MUSIC spectrum in (23), both the eigenvectors e′u and the search vector a′nm(Φ ) are real, therefore, the search function is real and the method is called real SHMUSIC. The real computations can save a considerable amount of computations. According to the following phase shift characteristic [15] of the spherical harmonic

Y nm(π − ϑ, π + φ) = ( − 1)nY nm(ϑ, φ ).

4

(23)

where a′nm(Φ ) = Q¯ a¯ nm(Φ ) with the real element

⎧ 2 Re{a¯ m(Φ)} , m < 0, n ⎪ ⎪ m a n′ (Φ) = ⎨ a¯ nm(Φ), m = 0, ⎪ m ⎪ ¯ Φ − { ( )} > 0, Im a m 2 , ⎩ n

SHMUSIC RVMUSIC LC-MUSIC,g=13 CRB

(25)

The elements of the real search vector have the following property

⎧ −a¯ ′m (ϑ, φ), n a¯ ′m n (ϑ , π + φ ) = ⎨ ⎪ m ⎩ a¯ ′n (ϑ, φ), ⎪

m

is

m

odd,

is

even,

(26)

and

⎧ −a¯ ′m (ϑ, φ), n a¯ ′m n (π − ϑ , φ ) = ⎨ ⎪ m ⎩ a¯ ′n (ϑ, φ), ⎪

n+m n+m

is is

odd, even.

(27)

According to (26), we can select the observation SFT coefficients corresponding to the search vector where m is even. Consequently when performing the spectrum search, we can just search the azimuth area of 0 ≤ φ ≤ π . Then we compute P (ϑ, φ) and P (ϑ, π + φ). If there is only one source, P (ϑ, φ) or P (ϑ, π + φ) must be much larger than the other. The larger value corresponds to the true DOA. If there are two approximately symmetrical sources,

Q. Huang et al. / Signal Processing 131 (2017) 441–446

1.8 SHMUSIC RVMUSIC LC-MUSIC,g=13 CRB

RM SE of elev ation (degree)

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

SNR (dB) Fig. 6. RMSE of elevation estimation versus SNR. Table 3 The spectrum results of P(39°, 40°) and P(141°, 40°) for one source and two sources. Simulation

P (39°, 40°)

P (141°, 40°)

One Source Two Sources

0.994 0.935

0.005 0.864

both P (ϑ, φ) and P (ϑ, π + φ) are large values (usually larger than 0.5). Based on (27), we can select the observation SFT coefficients corresponding to the search vector where m þn is even. When performing the spectrum search, we can just search the elevation area of 0 ≤ ϑ ≤ π /2. Then we compute P (ϑ, φ) and P (π − ϑ, φ) to determine one source or two sources. Considering both (26) and (27), we can select the observation SFT coefficients corresponding to the search vector where both m and m þn are even. In this situation, we can just search the elevation area of 0 ≤ ϑ ≤ π /2 and the azimuth area of 0 ≤ φ ≤ π . Then we compute P (ϑ, φ), P (π − ϑ, φ), P (ϑ, π + φ), and P (π − ϑ, π + φ) to determine the number of sources and their DOAs. These low complexity RSHMUSIC (LC-RSHMUSIC) methods further reduce the computation cost largely. We use g to denote the number of the observation SFT coefficients dependent on the values of n and m. As the three LC-RSHMUSIC methods choose different g, they have different DOA estimation performance.

4. Simulations In this section, some simulations are conducted in free-field environment to demonstrate the performance of the proposed methods. A 32-element uniform sampling of spherical array with the radius R¼4.2 cm is chosen to estimate source locations. The maximum order of the spherical harmonics is 4. In order to show the resolution power of our proposed methods, we compare them with the SHMUSIC method in (7). As the maximum order of the spherical harmonics is 4, the number g of the three LC-RSHMUSIC methods is 13, 15, and 9, respectively. Table 2 shows the relationship between n and m of the observation SFT coefficients and the steering vectors of the three LC-RSHMUSIC methods. Here we compare the LC-RSHMUSIC when g is 13, RSHMUSIC, and SHMUSIC methods. The SNR is set as 5 dB and the

445

number of snapshots is set as 200. Assuming two uncorrelated farfield narrowband signals with the same elevation impinge on the spherical array from (65°, 46°) and (65°, 59°), the MUSIC spectrum versus azimuth is shown in Fig. 1. Then the azimuths of two signals are kept the same and elevations are different. Two uncorrelated far-field narrowband signals impinge on the spherical array from (39°, 40°) and (51°, 40°). The result is shown in Fig. 2. These results indicate that the RSHMUSIC and LC-RSHMUSIC ap^ proaches have a better resolution than the SHMUSIC method. R nm is a forward only averaged estimate for the covariance matrix. According to [2], the matrix R˜ obtained by block forward-backward (FB) averaging is the optimal block centrohemitian estimator ^ of R nm in the sense of Euclidean distance. Therefore, the RSHMUSIC and LC-RSHMUSIC methods have the improved performance both in the estimation accuracy and resolution. The probability of correct estimation is used to evaluate the DOA estimation performance with 200 independent Monte Carlo trials. We compare the LC-RSHMUSIC when g is 13, 15, and 9, RSHMUSIC, and SHMUSIC methods. The number of snapshots is fixed at 200 and SNR changes from  5 dB to 5 dB. Two uncorrelated far-field narrowband signals impinge on the spherical array from (55°, 39°) and (55°, 58°), respectively. The angle search step is 0.1°. We consider the estimation right if max d = 1,2{|φ^d − φd|} is less than 1°, where φ^d stands for the estimated azimuth for the dth signal. The probability of correct estimation for azimuth is shown in Fig. 3. Assuming two signals impinge on the spherical array from (30°, 40°) and (52°, 40°), the probability of correct estimation for elevation is shown in Fig. 4. In this condition, the two ^ − ϑ |} is less than 1°, signals are estimated correctly if max {|ϑ d = 1,2

d

d

^ stands for the estimated elevation for the dth signal. where ϑ d From Figs. 3 and 4 we can know that the RSHMUSIC method has a better performance than the SHMUSIC method. The LC-RSHMUSIC methods are a little worse than the RSHMUSIC, but are better than SHMUSIC. Root mean square error (RMSE) is adopted to evaluate the DOA estimation performance. Cramer-Rao Bound (CRB) was analyzed for DOA estimation in spherical harmonics domain [19]. It provides a lower bound on the variance of estimated DOA and defines the ultimate accuracy. Two uncorrelated far-field narrowband signals impinge on the spherical array from (35°, 35°) and (35°, 57°). The number of snapshot is set as 200 with SNR varying from 0 dB to 20 dB. The RMSEs of 200 independent Monte Carlo trials are shown in Fig. 5. From Fig. 5 we can see that the RMSE of azimuth using the low complexity method is a little larger than that using RSHMUSIC method. But it is better than SHMUSIC method when SNR is low. Then we assume two signals impinge on the spherical array from (30°, 40°) and (52°, 40°). The RMSEs of elevation are shown in Fig. 6. These results indicate that the LC-RSHMUSIC method has a better DOA estimation performance than SHMUSIC when SNR is low. In order to illustrate the low computation amount of our methods, we show the CPU time of these methods when using the computer with Intel i7 CPU, 4 G RAM. The EVD of SHMUSIC, RSHMUSIC and LC-RSHMUSIC when g ¼13 costs 0.49 ms, 0.27 ms, and 0.1 ms, respectively. The time of one calculation of spectrum is 0.14 ms, 0.067 ms, and 0.052 ms, respectively. From these results we can know that our methods have a lower computation amount. In order to illustrate the DOA estimation for the signals in the symmetrical location clearly, we conducted two simulations by assuming one signal impinges on the array from (39°,40°) and two signals from (39°, 40°) and (141°, 40°). The SNR and the number of snapshots are set as 0 dB and 200, respectively. The spectrum values of P(39°, 40°) and P(141°, 40°) are shown in Table 3. When there is only one source, the spectrum value in the true DOA is much larger than that in its symmetrical mirror. While there are

446

Q. Huang et al. / Signal Processing 131 (2017) 441–446

two sources in the symmetrical locations, both of the two spectrum values are large. Therefore, LC-RSHMUSIC can give the correct number of sources and estimate true DOAs.

5. Conclusions In this paper, we propose two unitary transformations for DOA estimation using spherical arrays. Based on the unitary transformations, a real SHMUSIC method with lower amount of computations is proposed to estimate the signals locations. The performance is improved when the SNR is low because the new covariance matrix exploits the block FB averaging and block centrohermitian property. Moreover, the computation cost can be further reduced by exploiting the phase shift characteristic of the real search vector. This method can be extended in many direction finding and spectral estimation problems.

Acknowledgment The authors would like to thank the editor and anonymous reviewers for their valuable comments. This work was supported by National Natural Science Foundation (61571279, 61501288) and Shanghai Natural Science Foundation (14ZR1415000) of China.

Appendix Using the block centrohermitian property C = ¯JC*¯J , ¯¯ JJ = IM , and ¯ * Ξ J = Ξ, the conjugate of ΞCΞH can be expressed as

(ΞCΞH )* = Ξ*C*ΞT = Ξ*¯¯ JJC*¯¯ JJ ΞT = ΞCΞH .

(28)

H

Therefore, ΞCΞ is real.

References [1] R.O. Schmidt, Multiple emitter location and signal parameter estimation, IEEE

Trans. Antennas Propag. 34 (3) (1986) 276–280. [2] K.C. Huarng, C.C. Yeh, A unitary transformation method for angle-of-arrival estimation, IEEE Trans. Signal Process. 39 (1991) 975–977. [3] M. Pesavento, A.B. Gershman, M. Haardt, Unitary root-MUSIC with a real-valued eigendecomposition: a theoretical and experimental performance study, IEEE Trans. Signal Process. 48 (5) (2000) 1306–1314. [4] N. Yilmazer, T.K. Sarkar, 2-D unitary matrix pencil method for efficient direction of arrival estimation, Digit. Signal Process. 16 (6) (2006) 767–781. [5] J. Jiang, L. Gan, Decoupled unitary esprit algorithm for 2-D DOA estimation, Prog. Electromagn. Res. 29 (2012) 219–234. [6] I. Kazemi, M.R. Moniri, R.S. Kandovan, Optimization of angle-of-Arrival estimation via real-valued sparse representation with circular array radar, IEEE Access 1 (2013) 404–407. [7] X. Li, S. Yan, X. Ma, C. Hou, Spherical harmonics MUSIC versus conventional MUSIC, Appl. Acoust. 72 (2011) 646–652. [8] H. Sun, E. Mabande, K. Kowalczyk, W. Kellermann, Localization of distinct reflections in rooms using spherical microphone array eigenbeam processing, J. Acoust. Soc. Am. 131 (4) (2012) 2828–2840. [9] Q. Huang, T. Wang, Acoustic source localization in mixed field using spherical microphone arrays, EURASIP J. Adv. Signal Process 90 (2014). [10] L. Kumar, B. Guoan, R.M. Hegde, The spherical harmonics root-MUSIC, in: Proceedings of 41st IEEE International Conference on Acoustic Speech and Signal Processing (ICASSP), Shanghai, China, 2016, pp. 3046–3050. [11] F. Yan, M. Jin, S. Liu, X. Qiao, Real-Valued MUSIC for efficient direction estimation with arbitrary array geometries, IEEE Trans. Signal Process. 62 (6) (2014) 1548–1560. [12] R. Goossens, H. Rogier, Unitary spherical ESPRIT: 2-D angle estimation with spherical arrays for scalar fields, IET Signal Process. 3 (3) (2009) 221–231. [13] S. Yan, H. Sun, X. Ma, U.P. Svensson, C. Hou, Time-domain implementation of broadband beamformer in spherical harmonics domain, IEEE Trans. Audio Speech Lang. Process. 19 (5) (2011) 1221–1230. [14] B. Rafaely, Analysis and design of spherical microphone arrays, IEEE Trans. Speech Audio Process. 13 (1) (2005) 135–143. [15] E.G. Williams, Spherical Waves Fourier Acoustics, Chapter 6, Academic press, New York 1999, pp. 183–197. [16] S. Yan, H. Sun, U.P. Svensson, X. Ma, J.M. Hovem, Optimal modal beamforming for spherical microphone arrays, IEEE Trans. Audio Speech Lang. Process. 19 (2) (2011) 361–371. [17] A. Lee, Centrohermitian and skew-centrohermitian Matrices, Linear Algebra Its Appl 29 (1980) 205–210. [18] D.A. Linebarger, R.D. DeGroat, E.M. Dowling, Efficient direction-finding methods employing forward/backward averaging, IEEE Trans. Signal Process. 42 (8) (1994) 2136–2145. [19] L. Kumar, R.M. Hegde, Stochastic Cramer-Rao bound analysis for DOA estimation in spherical harmonics domain, IEEE Signal Process. Lett. 22 (8) (2015) 1030–1034.