社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

高次 Root-Music とその実多項式形
曹祥†

上手 洋子†

西尾 芳文†

辛景民††† 徳島大学工学部 〒 770–8506 徳島県徳島市南常三島町 2-1†† Xi’an Jiaotong University No.28 Xianning West Road, Xi’an 710049, China
E-mail: †{caoxiang, uwate, nishio}@ee.tokushima-u.ac.jp, †† jxin@mail.xjtu.edu.cn
あらまし

本研究では，Root-MUSIC 到来方向推定法の実多項式を提案する．共形変換を用いて Root-MUSIC を実係数多項式への変換を行う．そして，上半単位円上の変数を用いた高次 Root-MUSIC の提案を行う．提案手法は J. Selva

2005 の手法より優れた性能を示した．
キーワード パラメータ概算，ULA，アレイ信号処理

A higher order Root-MUSIC and its real polynomial form
Xiang CAO† , Yoko UWATE† , Yoshifumi NISHIO† , and Jingmin XIN††† Department of Electrical and Electronic Engineering, Tokushima University
2-1 Minami-Josanjima, Tokushimashi, Tokushima, 770-8506, Japan†† Department of Institute of Artificial Intelligence and Robotics,
Xi’an Jiaotong University No.28 Xianning West Road, Xi’an 710049, China
E-mail: †{caoxiang, uwate, nishio}@ee.tokushima-u.ac.jp, †† jxin@mail.xjtu.edu.cnAbstract This paper presents a real polynomial form of the popular Root-MUSIC direction-of-arrival (DOA)
algorithm. To transform Root-MUSIC into real coefficients polynomial with the help of conformal transformation,
we propose a higher order Root-MUSIC whose variable is only defined on the unit upper semicircle. Furthermore,this paper also improves the results obtained by J.Selva in 2005.
Key words Parameter estimation, uniform linear array (ULA), array signal processing

1. Introduction

mial in [4] has the same degree with standard Root-MUSIC.
In this paper, we improve the result and give the real poly-

The multiple signal characterization (MUSIC) algorithm

nomial form of Root-MUSIC in [−π, π]. The same mapping

for direction finding [1] has been advanced for many years.

idea also appears in [5] to estimate single target. In [5], theTo a certain degree, the more search points one use, theinterested range is mapped into real line by a tangent func-

more accurate results one will obtain. The main computa-tion, while the transformation function is not monotonic overtional burden locate the process of spectral search. Thus,target range. That is to say, the target may not be uniquely

many efficient methods are proposed to reduce the search

determined from the inverse transformation.

complexity of spectral MUSIC.

In this paper, we propose a higher order Root-MUSICTo avoid the tremendous computation load of spectral

whose variable is only defined on the unit upper semicir-

search, Root-MUSIC [2], [3] was proposed. The search-freecle. The similar Root-MUSIC then is transformed into real

method has an improved resolution threshold relative to con-polynomial form. Our method treats the whole range of in-ventional MUSIC. Taking advantage of conformal transfor-terest in uniform linear array (ULA). The process of solving

mation, the author of [4] translated the Root-MUSIC into an

this polynomial just involves real arithmetic [6]. This method

univariate polynomial with real coefficients. Here, we havealso can be used in unitary Root-MUSIC [7] and unitary MU-

to point out that the interested spatial frequency in Root-

SIC [8].

MUSIC lies in the interval [−π, π] instead of [0, π]. Since only
half of field-of-view is considered, the real coefficients polyno—1—

is the number of spectral search points. Nevertheless, the2. Data Model

subspace decomposition in (5) costs O(M 2 p) complex flops

Assume p far-field narrowband signals {sk } impinging on[11], [12]. We can conclude that the main computational loada ULA with M (M > p) sensors. The output vector y(t) ofof spectral MUSIC is the search steps. Consequently, it is

the array at time t can be written asnecessary to reduce the complexity of spectral search process.

y(t) = A(θ)s(t) + n(t),(1)Next, we will deduce the real polynomial form of RootMUSIC with the help of conformal transformation. By defin-where s(t) is the vector of incident signals and n(t) is theing a translation

vector of additive noises. The so-called array steering matrixz = ej( 2 + 2 ) ,and steering vector have the following form:
A(θ) = [a(θ1 ), a(θ2 ), · · · , a(θp )22]hiT
2π2πa(θi ) = 1, ej λ d sin(θi ) , · · · , ej λ (M −1)d sin(θi ) ,here, (·)

T(2)

ωπ(8)

the real point ω belonging to the domain [−π, π] can bemapped to the unique complex point z belonging to the unit
upper semicircle. Furthermore, let us introduce a conformal(3)

denotes transpose operator, d and λ stand forthe distance between two consecutive sensors and the iden-transformationu = −jz−j.
z+j

(9)

tical wavelength for all signals, respectively. The restriction

The above function is a one-to-one mapping that takes theof Directions-of-arrival (DOAs) θ1 , · · · , θp lie in the interval − π2 , π2 for the purpose of eliminating the ambiguity ofunit upper semicircle to the interval [−1, 1] of real line [4],ULA [9].

Next, we define spatial frequency ωi = 2πd sin(θi ).
λ(3) becomeshiTa(ωi ) = 1, ejωi , · · · , ej(M −1)ωi .

Then,(4)[5], [13]. Substituting (8) into (9), we get the relationshipbetween ω ∈ [−π, π] and u ∈ [−1, 1] as follows:ω.
u = tan4(10)This is a monotonic function in [−π, π] and its range is[−1, 1]. Thus, once the variable u is estimated, the unique

For guaranteeing the uniqueness of a(ωi ), i = 1, 2, · · · , p, thespatial frequency ω can be obtained by ω = 4 arctan u. For

spatial frequency ωi are in [−π, π] instead of [0, π] in [4]. In

this purpose, the array vector response a(ω) must be con-other words, the author in [4] only consider the half of thenected to u.

DOA range of interest. Since the one-to-one correspondence
between ωi and θi , we only focus on the problem of estimating parameters ωi , i = 1, 2, · · · , p in this paper.

3. Real polynomial form of Root-MUSIC
Under the assumptions of the uncorrelated signals and theUsing z = −j u−jand ejω = −z 2 ( the two equations canu+j
be easily obtained from (8) and (9) ), we get[a(ω)]m=ej(m−1)ω = (−1)m−1 z 2(m−1)
2(m−1)u−j=, m = 1, 2, · · · , M.
u+j

(11)

white Gaussian noise, the array covariance matrix can be

Multiplying by (u + j)2(M −1) (here, we note u =| −j ), thedecomposed intoabove equation then becomesHR = E{y(t)y H (t)} = U s Λs U Hs + U n Λn U n .

(5)

It is clear that the columns of A(θ) have the maximum projection on the signal subspace U s and are orthogonal to the(u + j)2(M −1) [a(ω)]m=(u − j)2(m−1) (u + j)2(M −m)=cm,2M −2 u2M −2 + · · · + cm,0=cTm v(u), m = 1, 2, · · · , M(12)noise subspace U n . With this observation, two natural esti-

where the (2M − 1) × 1 coefficient vector cm = [cm,0 , cm,1 ,mation criterions are to find the local maxima or minima [1],· · · , cm,2M −2 ]T can be calculated by convolution operationviz.

[5], [14], v(u) = [1, u, · · · , u2M −2 ]T is unknown real vector.

ω̂ = arg max f1 (ω) = arg max aH (ω)U s U Hs a(ω)
ω
ω

(6)

Let us define an M ×(2M −1) matrix C = [c1 , c2 , · · · , cM ]T ,the array steering vector a(ω) can be expressedora(ω) = (u + j)2(1−M ) Cv(u).
ω̂ = arg min f2 (ω) = arg min aωωH(ω)U n U Hn a(ω).

(13)

(7)
In ULA, the Root-MUSIC always has more resolution com-

The total number of complex flops in (6) is Jp(2M +1) andpared with spectral MUSIC [3]. Thus, it is interesting to lookthe one in (7) is J(M −p)(2M +1) [10], where J (J ≫ M > p)for a real polynomial form of Root-MUSIC. Using z = ejω ,—2—ω ∈ [0, π] and (9), a real polynomial with degree 2(M − 1)0.655

has been obtained in [4]. However, If we consider the whole

0.65domain ω ∈ [−π, π], the conformal transformation equation

0.645(9) will not exist in z = ejω , but exists in (8).

0.64

Thus, we can define the mth element of vector a(z) as[a(z)]m = (−1)m−1 z 2(m−1) , m = 1, 2, · · · , M.

(14)

The function f2 (ω) in (7) then becomesf2 (z) = aH(z)U n U Hn a(z) = 0.

0.635
RMSE

follows:

Real polynomial
Real-Imag polynomial
Root-MUSIC

0.63

0.625

0.62(15)0.615This is a polynomial of degree 4(M − 1) with complex coeffi-0.61cients and the variable z is defined on the unit upper semicir-0.605-10-8-6-4-2

cle instead of the entire unit circle in the conventional RootMUSIC with 2(M − 1)-degree complex coefficients. Using
the equation a(z) = (u + j)
Tv (u)CH2(1−M )

2468

10

RMSE of spatial frequency ω=-0.3. The number of snapshots N = 100

Cv(u) , we get

U nU Hn Cv(u) = 0.

図 10SNR(16)The roots of polynomial in (15) appear in reciprocal pairs
z, z1∗ [9] and could not be on the unit circle because of the4. EXPERIMENTS
In this section, an examples was conducted by computer

effects of noise. The pair roots are related to u by conformal

simulation. In all scenarios, two equal-power uncorrelatedtransformation (9),signals are considered and a ULA have 9 sensors. Mean-z−j
1 − jz ∗u1 = −j, u2 = −j.
(17)z+j
1 + jz ∗Obviously, we have u1 = u∗2 . It is reasonable to conclude thatwhile, we set spatial frequency ω1 = −0.3, ω2 = 1.2. Here,the roots in (16) also appear in conjugate pairs like (u, u∗ ).
As a consequence, the equation (16) is a real coefficientspolynomial, namelynov T (u)Re C H U n U Hn C v(u) = 0.

ω1 can not be solved in [4].
Example: Performance of real polynomial based on RootMUSICIn this case, we mainly test the validity of real coefficientspolynomial (18). For this purpose, let us introduce the imag-(18)Next, we will analyze how to estimate the spatial frequency
ωi , i = 1, 2, · · · , p from the conjugate roots (ui , u∗i ) ,i = 1, 2, · · · , 2(M − 1) of the above polynomial. Assumingz = a + jb, b >= 0, the equation (9) becomes−2a − j a2 + b2 − 1
a + jb − j
u = −j=.
(19)
a + jb + j
a2 + (b + 1)2

inary part of matrix H = C H U n U Hn C. After simple derivation, we obtainv T (u) {Re {H} + Im {H}} v(u) = 0.

(20)This approach using (20) named real-imag polynomial andwas indicated by red line in Figure. 1. It is worth pointing
out that (18) and (20) give the same results. That is to say,When the roots z in (15) are close to the unit circle, thethe coefficients of polynomial (16) are real. Meanwhile, theimaginary part of u must tend to zero. Therefore, p real

behavior of our method have been improved a little at lower

parts of ui can be picked up through comparing the absolu-

signal-to-noise ratio (SNR) than Root-MUSIC [2].

tion value of imaginary part of ui , i = 1, 2, · · · , 2(M −1). The

5. CONCLUSIONSspatial frequency will be given by ωi = 4 arctan [Re (ui )] , i =For improving the results in [4], we propose a higher order

1, 2, · · · , p.
Note that, although the polynomial in (18) have 4(M − 1)-

Root-MUSIC whose variable is only defined on the unit up-degree, its coefficients are real and roots appear in conjugate

per semicircle. Then, the complex equation is transformedpairs. Thus, we only need to solve 2(M − 1) roots. For real

into real polynomial form with the help of conformal trans-polynomial, the process of solution can avoid all complex

formation. The process of solution can avoid all complex

arithmetic using Bairstow’s method [6]. In practice, finding

arithmetic.

roots of a function usually involve iteration. The operation
2in (18) is O(16M ) real flops per iteration, while the RootMUSIC requires O(24M 2 ) real flops per iteration (1 complex flop = 6 real flops) even though polynomial have degree2(M − 1).

References[1] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag., vol.
53, pp. 276–280, Mar. 1986.
[2] A. J. Barabell, “Improving the resolution performance ofeigenstructure-based direction-finding algorithms,” in Proc.

—3—

ICASSP, Boston, MA, May 1983, pp. 336–339.
[3] B. D. Rao and K. V. S. Hari, “Performance analysis of rootmusic,” IEEE Trans. Acoust.,Speech, Signal Process., vol.
37, pp. 1939–1949, Dec. 1989.
[4] J. Selva, “Computation of spectral and root music throughreal polynomial rooting,” IEEE Trans. Signal Process., vol.
53, pp. 1923–1927, May 2005.
[5] J. X. Wu, T. Wang, and Z. Bao, “Fast realization of maximum likelihood angle estimation with small adaptive uniform linear array,” IEEE Trans. Antennas Propag., vol. 58,pp. 3951–3960, Dec. 2010.
[6] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.
Flannery, Numerical recipes in C, Cambridge Univ. Press,Cambridge, U.K., 1997.
[7] M. Pesavento, A. B. Gershman, and M. Haardt, “Unitaryroot-music with a real-valued eigendecomposition: a theoretical and experimental performance study,” IEEE Trans.
Signal Process., vol. 48, pp. 1306–1314, May 2000.
[8] K. C. Huarng and C. C. Yeh, “A unitary transformation
method for angle-of-arrival estimation,” IEEE Trans. signal
process., vol. 39, pp. 975–977, Apr. 1991.
[9] P. Stoica and R. L. Moses, Spectral analysis of signals,
Prentice-Hall, 2005.
[10] D. S. Watkins, Fundamentals of matrix computations, Wiley, 2004.
[11] Y. Zhang and B. P. Ng, “Music-like doa estimation without estimating the number of sources,” IEEE Trans. Signal
Process., vol. 58, pp. 1668–1676, Mar. 2010.
[12] F. G. Yan, M. Jin, and X. Qao, “Low-complexity doa estimation based on compressed music and its performanceanalysis,” IEEE Trans. Signal Process., vol. 61, pp. 1915–1930, Apr. 2013.
[13] L. Blanco, J. Serra, and M. Nájar, “Conformal transformation for efficient root-toa estimation,” in Proc. ISCCSP,Malta, Mar. 2008, pp. 1314–1319.
[14] M. Costa and V. Koivunen, “Incorporating array nonidealities into adaptive capon beamformer for source tracking,”
in Proc. SAM, Hoboken,NJ, Jun. 2012, pp. 445–448.

—4—