International Workshop on Acoustic Signal Enhancement 2012, 4-6 September 2012, Aachen

UNDERDETERMINED DOA ESTIMATION BY THE NON-LINEAR MUSICEXPLOITING HIGHER-ORDER MOMENTS
Yuya Sugimoto1 , Shigeki Miyabe1 , Takeshi Yamada1 , Shoji Makino1 and Biing-Hwang (Fred) Juang2
1

2

University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki, 305-8577 Japan
Georgia Institute of Technology, 777 Atlantic Dr. NW. Atlanta, GA 30332-0250 USA
Email: sugimoto@mmlab.cs.tsukuba.ac.jp, {miyabe, makino}@tara.tsukuba.ac.jp,
takeshi@cs.tsukuba.ac.jp, juang@ece.gatech.edu
ABSTRACTThis paper describes a new approach for extending MUltiple SIgnal Classification (MUSIC) to underdetermined direction-of-arrival(DOA) estimation with high resolution by exploiting higher-ordermoments. The proposed method maps the observed signals nonlinearly onto a space of expanded dimensions, in which signal statistics
are analyzed. The covariance matrix in the higher-dimensional spacecorresponds to the higher-order cross moment matrix in the original
space of the observed signals. Since the dimensionality of the noise
subspace is increased by the mapping, the proposed method achieves
higher resolution DOA estimation than the standard MUSIC, andalso offers the ability to estimate DOAs in underdetermined conditions. We compared the characteristic of the proposed method withthat of the conventional 2q-MUSIC utilizing higher-order cumulants
theoretically and experimentally.
Index Terms— MUSIC, Higher-order statistics, Microphone array, Underdetermined DOA estimation
1. INTRODUCTION
In array signal processing, direction of arrival (DOA) estimation
is essential in various applications such as source separation andnoise reduction. MUltiple SIgnal Classification (MUSIC) [1] is a
popular high-resolution DOA estimation method which employs thesubspace analysis of the observed signals. However, MUSIC is restricted by the dimensionality of the covariance matrix because MUSIC requires the identification of the noise subspace that is orthogonal to the transfer function vectors of the observed signals. To estimate N sound sources, M > N sensors are required and the DOA
estimation performance deteriorates as N approaches M . Thus,
MUSIC has a problem in that an increase in the array scale is unavoidable when estimating a large number of sources.
To overcome this issue, several extensions of MUSIC have been
proposed that increase the dimensionality of the observed signals
virtually by using higher-order cumulants. For example, a MUSIClike algorithm [2] exploits fourth-order cumulants, and its extension,
2q-MUSIC [3], exploits cumulants of an arbitrary even order. Thesemethods improve the resolution of DOA estimation and can also estimate DOAs in underdetermined conditions, where the number ofsources exceeds the number of sensors, by increasing the signal expressiveness by way of subspace expansion.
A speech enhancement method called complementary beamforming [4] has been proposed for scenarios with more noise sources
This work was supported by KAKENHI 23760325.

than sensors, and this approach was used for a DOA estimation problem [5]. Complementary beamforming has also been extended andexplained as the mapping of a signal onto higher dimensional space
using a kernel functions [6]. The approach used in these methods is
especially effective for MUSIC to which dimensionality is critical.
In this paper, utilizing mapping for MUSIC, we propose mapped
MUSIC for underdetermined DOA estimation. By increasing thedimensionality of the noise subspace with the higher-dimensional
mapping of the observed signal, the resolution of DOA estimation
is improved and DOA estimation in underdetermined conditions
is achieved. We describe a class of mapping suitable for mapped
MUSIC, which allows us to analyze the cross moments of arbitraryeven orders. Also we show an efficient way to calculate the crossmoments for fourth and sixth orders. Moreover, we discuss therelation between the cross moment used in the proposed mapped
MUSIC and the cross cumulant used in 2q-MUSIC, and show thatthe analysis of the cross moment is more computationally efficient.
Experimental results in a simulation of speech DOA estimation
reveals that the proposed mapped MUSIC can achieve high resolution DOA estimation similar to 2q-MUSIC but with much lesscomputational complexity.
2. PROBLEM STATEMENT
Throughout this paper, signals are expressed as complex amplitudes.
Observed signals can be modeled asx(ω, t) = [x1 (ω, t), · · · , xM (ω, t)]T=N∑ai (ω)si (ω, t) + n(ω, t),(1)n(ω, t) = [n1 (ω, t), · · · , nM (ω, t)]T ,(2)ai (ω) = [a1,i (ω), · · · , aM,i (ω)] ,

(3)

i=1

T

where t = 1, · · · , L is a time frame index, ω is the angular frequency, M is the number of sensors, N is the number of soundsources, [·]T denotes transpose, si (ω, t) is the complex amplitudeof the ith sound source, xj (ω, t) is the complex amplitude of thesignal observed with the jth sensor, nj (ω, t) is the complex amplitude of the noise observed with the jth sensor, and aj,i (ω) denotes
the transfer function from the ith source to the jth sensor.
In the signal model expressed by Eq. (1), each sensor observes a mixture of source signals and noise signals. The problem in this paper is to estimate the DOAs of the source signals s1 (ω, t), · · · , sN (ω, t) by finding steering vectorsb(ω; θ1 ), · · · , b(ω; θN ) similar to the transfer function vectors
a1 (ω), · · · , aN (ω).

3. PROPOSED METHOD
3.1. DOA estimation algorithm of mapped MUSIC
Mapped MUSIC maps the M -dimensional observed vector x(ω, t)onto M 0 -dimensional Euclidian space (M 0 ≥ M ) with a nonlinear0function φ : CM → CM and conducts a similar analysis to MUSIC
with vectors φ(x(ω, t)). To estimate DOAs properly with mapped
MUSIC, the information about the correlations between each microphone signal must be retained after mapping. This requirement isexpressed as the following three conditions:||x||2 ≥ ||y||2 → ||φ(x)||2 ≥ ||φ(y)||2
2. The origin remains intact.

(4)x = 0 → φ(x) = 0
3. The orthogonality between vectors is preserved.

(5)

xH y = 0 → φH (x)φ(y) = 0

R(ω) = E[φ(x(ω, t))φH (x(ω, t))],(7)

where E[·] denotes the expectation of the argument, and [·]H denotes
a complex conjugate transpose. The following equations are givenby the eigen decomposition of covariance matrix R(ω),R(ω) = V(ω)E(ω)VH (ω),(8)

V(ω) = [v1 (ω), · · · , vM 0 (ω)], V (ω)V(ω) = IM 0 , (9)E(ω) = diag[e1 (ω), · · · , eM 0 (ω)],(10)
e1 (ω) ≥ · · · ≥ eM 0 (ω),
H

M 0 = dim[φ(x(ω, t))],(11)

where v1 (ω), · · · , vM 0 (ω) are the eigenvectors for each eigenvalue
e1 (ω), · · · , eM 0 (ω), respectively, Ii denotes the i-dimensional identity matrix, diag[·] is a diagonal matrix consisting of elements within
the argument vector, and dim[·] is the dimensionality of the argument. With N 0 denoting the number of large eigenvalues corresponding to the mappings of the directional sound source signals,
we define the spanned subspace by v1 (ω), · · · , vN 0 (ω) as the signal subspace S(ω) in the mapping space. Incidentally, this paper
assumes N 0 = N .
(12)where span[·] denotes subspace spanned by argument vectors.
Moreover, its orthogonal complement in span[v1 (ω), · · · , vM 0 (ω)]
is defined as the noise subspace. Now, the following relation exists
between the transfer function vectors a1 (ω), · · · , aN (ω) and thevectors vN 0 +1 (ω), · · · , vM 0 (ω) in the noise subspace.
φH (ai (ω))vj (ω) = 0,(i, j) ∈ {i = 1, · · · , N j = N 0 + 1, · · · , M 0 }.

(13)

Let us define the following MUSIC score of mapped MUSIC f (ω; θ)utilizing the orthogonality in Eq. (13).
1

φd (x) =

[ d
∏

xl⊕
c1l , · · · ,

l=1

,

(14)

d
∏

]T
xl⊕
cM d l

,

(16)

l=1

x = [x1 , · · · , xM ]T ,{ ∗x
(if l is odd)
xl⊕ =
,
x
(if l is even)

(17)
(18)

where {·}∗ denotes a complex conjugate. Then ckl is an index specifying the element number of vector x for calculating φd , and thisis the (k, l) element of the M d × d index matrix C that arranges drepeated-permutations of M in an arbitrary order as its row vectors:
1 1 ··· 1
 1 1 ··· 2 
C = [ckl ]kl = 
(19)
,
..


.
M M ··· M
where [·]ij denotes a matrix consisting of the argument as its (i, j)element. In practice, the covariance matrix of the mapping of the observed signals φd (x(ω, t)) is expressed as a 2d-th order cross moment matrix according to the following equation,
Rd (ω) = E[φd (x(ω, t))φHd (x(ω, t))] = [rij ]ij ,(i, j = 1, · · · , M d ),
)∗ ]
[( d
)( d
∏ ⊕
∏ ⊕
rij = E
xcil
xcjl
.
l=1

S(ω) , span[v1 (ω), · · · , vN 0 (ω)],|φH (b(ω; θ))vi (ω)|(15)

As shown in Sect. 3.1, mapped MUSIC can use any mapping function satisfying Eqs. (4–6), but its property is altered by the choice of
mapping. In this paper, to evaluate the properties of mapping quantad
tively, we focus on the following mapping φd : CM → CM , which
gives a 2d-th order cross moment matrix as its covariance matrix.

(6)

We describe the DOA estimation algorithm of mapped MUSIC
with the mapping φ satisfying Eqs. (4–6). The covariance matrix ofφ(x(ω, t)) is expressed asi=N 0 +1

φ1 (x) = x.
3.2. Mapping for analysis of 2d-th order moments

1. The magnitude relation of the norm is retained.

f (ω; θ) , ∑M 0

where b(ω; θ) (||b(ω; θ)||2 = 1) is a steering vector for direction θ. The score f (ω; θ) takes a high value when the θvalue is close to φH (b(ω; θ))φ(ai (ω)) ≈ ||φ(ai (ω))||2 forany i = 1, · · · , N . Then, the conditions (4–6) result in relationbH (ω; θ)ai (ω) ≈ ||ai (ω)||2 . Mapped MUSIC utilizes this property
and estimates DOA θ to make the value of f (ω; θ) high. Here, notethat mapped MUSIC with the following mapping with no signal
modification is equivalent to MUSIC.

(20)
(21)

l=1

We can see the increase in the dimensionality of the covariance matrix from M to M d by the definition of the mapping φd : CM →d
CM . Then mapped MUSIC enables us to estimate the DOAs in underdetermined conditions due to the enhanced expressiveness of thenoise subspace derived from the increase in dimensionality.
Equation (16) gives the mapping φd for the analysis of the evenorder cross moment, but these mappings can be replaced by an arbitrary mapping φ0d which gives the equivalent inner product:0H0φH
d (x)φd (y) = φd (x)φd (y),(22)

and there are more useful mappings than that given by Eq. (16) fromvarious viewpoints. For example, when the degree of the mapping2function is two, the following mapping φ02 : CM → RM satisfying Eq. (22) simplifies the calculation of the covariance matrix

and eigenvalue problem because its covariance matrix becomes realvalued.

have sharper peaks in the MUSIC score than mapped MUSIC if there
are a sufficient number of observed signals. However, as we show in
our experiments, mapped MUSIC performs as well as 2q-MUSIC inpractical conditions with a realistic number of observed signals.

φ02 (x) , [φ0abs (x)T , φ0re (x)T , φ0im (x)T ]T ,(23)

φ0abs (x) , [∀|xi |2 |1 ≤ i ≤ M ]T ,(24)

φ0re (x) ,2[∀Re[xi x∗j ]|2 ≤ i ≤ M, 1 ≤ j ≤ i − 1]T ,(25)

4. EXPERIMENT

2[∀Im[xi x∗j ]|2 ≤ i ≤ M, 1 ≤ j ≤ i − 1]T .

(26)

This section describes a simulation experiment to verify the efficiency of the proposed mapped MUSIC. We implemented MUSICand 2q-MUSIC (q = 2, 3) for comparison.

φ0im (x) ,√√When the degree is three, there are several dependent elements withthe same values in φ3 so that the rank of R03 becomes the same or32less than M +M
. Thus we can use the following mapping φ03 :2
M 3 +M 2

which satisfies Eq. (22) and collects the same
CM → C 2
elements. Such mapping simplifies the calculation of the covariancematrix and the eigenvalue problem by omitting dependent rows andcolumns from the covariance matrix.
φ03 (x) , [φ03a (x)T , φ03b (x)T , φ03c (x)T , φ03d (x)T ]T ,(27)

φ03a (x) , [∀|xi |2 x∗i |1 ≤ i ≤ M ]T ,(28)

T
φ03b (x) , [∀xi x∗2
j |1 ≤ i ≤ M, j 6= i] ,√φ03c (x) , 2[∀|xi |2 x∗j |1 ≤ i ≤ M, j 6= i]T ,√φ03d (x) , 2[∀xi x∗j x∗k |1 ≤ i ≤ M, j 6= k 6= i]T .

(29)
(30)

3.3. Comparison with 2q-MUSIC exploiting cumulants
This section discusses the differences between mapped MUSIC
based on an analysis of even order cross moments and 2q-MUSIC
based on an analysis of even order cross cumulants. The calculationof ith-order cumulants requires diverse multiplications of the moments, whose order is the same or less than i, and the computational
cost will increase rapidly as the order of analysis increases. Thus
mapped MUSIC, which can obtain a cross moment matrix with asingle mapping, is simpler than 2q-MUSIC and can more easily
expand the statistics of analysis to a higher order. As an example,
with the straightforward analysis of fourth and sixth order statistics,
these methods require the computational costs shown in Table 1,
where M is the number of sensors, L is the number of time framesof observed signals, and Dim. denotes dimension.
Table 1. Computational costs
Multiplication (times)
Eigen decomposition

2q-MUSIC (q = 3)
mapped MUSIC(d = 2)
mapped MUSIC(d = 3)

We conducted an evaluation of the DOA estimation and computational complexity of each method by using a simulated recording ofa speech mixture with a circular microphone array. We establishedthree different reverberant conditions by using the image method [7]with different reflection coefficients. To emulate noisy observation,diffused pink noises [8] with three different SNRs were added to theobserved signals, and to determine the influence of the reflection andnumber of sample frames, we also evaluated the performance withtwo different lengths of observed signals. Table 2 shows the otherexperimental conditions.

(31)

We evaluate mapped MUSIC with mappings φ02 and φ03 in Sect. 4.

2q-MUSIC (q = 2)

4.1. Experimental condition

L(18M 4 − 16M 2 )complex-valued matrix of M 2 Dim.
L(210M 6 − 206M 3 )complex-valued matrix of M 3 Dim.
L
(M 4 + 5M 2 − 2M )2real-valued matrix of M 2 Dim.
L
(M 6 + 2M 5 + M 4 + 34M 3 − 6M 2 )83
2complex-valued matrix of M +M
Dim.
2

While the 2q-th order cross cumulant matrix in 2q-MUSIC usually has the full rank of M q , the dimensionality of the mapping φdwhen d > 2 is less than M d as discussed in Sect. 3.2. In addition,
there is often a further reduction in the rank of the cross momentmatrix because of the environment. As a result, 2q-MUSIC tends toTable 2. Experimental conditions
Sensor array form
Sound sources# of sound sources# of sensors
Room size
Reverberation
Duration
Frame length
Frame shift length
Window function

3, 5
4Circular array with radius of 0.1 m
Speakers 1.5 m apart from array
Sampling frequency
16 kHzSNR
0,10,20 [dB]
5m×5m
T60 of 0,0.12,0.3 [s]
3,5 [s]1024 samples
512 samples
Hamming window

For the evaluation, we employed the root mean squared error(RMSE) between the estimated direction and the true sound direction. We evaluated 100 combinations of the positions of three or fivesources selected randomly from the directions {0◦ , 30◦ · · · , 330◦ }.
When the MUSIC score had fewer peaks than sources, we added apenalty equal to the average error for all directions. And to compare
the computational complexity, we recorded processing time of eachmethod in every trial.
DOA estimation with MUSIC when M ≤ N is performed byregarding the one-dimensional subspace concurrent with the minimum eigenvalue as the noise subspace in every frequency bin.
4.2. Experimental results
Figures 1–4 show the experimental results obtained under different conditions. The results reveal that the proposed mapped MUSIC and the 2q-MUSIC perform similarly better than MUSIC under all conditions and demonstrate the high performance of DOA
estimation even in underdetermined conditions. Under most conditions, mapped MUSIC with d = 2 performs as well as 2q-MUSIC.
Mapped MUSIC with d = 3 shows excellent accuracy in terms of
DOA estimation with a high SNR and is the method that performs
the best in this experiment. On the other hand, 2q-MUSIC showsslightly poorer performance, although 2q-MUSIC can be expected

102

102

MUSIC

MUSIC
2q-MUSIC (q=2)
2q-MUSIC (q=3)

100Mapping MUSIC(d=2)

RMSE [degree]

RMSE [degree]101
2q-MUSIC (q=2)

101

2q-MUSIC (q=3)
Mapping MUSIC
(d=2)

100

10-1
Mapping MUSIC(d=3)

Mapping MUSIC
(d=3)

10-2
SNR 0dB 10dB 20dB

T60=0 sSNR 0dB 10dB 20dB

10-1
SNR 0dB 10dB 20dB

SNR 0dB 10dB 20dB

T60=0.12 s

T60=0 s

T60=0.3 s

Fig. 1. Experimental results for 3 sources within 3 seconds.

SNR 0dB 10dB 20dB

SNR 0dB 10dB 20dB

T60=0.12 s

T60=0.3 s

Fig. 3. Experimental results for 5 sources within 3 seconds.

102

102
MUSIC

MUSIC

2q-MUSIC (q=3)

100Mapping MUSIC
(d=2)

10-1

RMSE [degree]

2q-MUSIC (q=2)

2q-MUSIC (q=2)

101

2q-MUSIC (q=3)
Mapping MUSIC
(d=2)

100

Mapping MUSIC
(d=3)

10-2
SNR 0dB 10dB 20dB

T60=0 sSNR 0dB 10dB 20dB

SNR 0dB 10dB 20dB

T60=0.12 s

Mapping MUSIC
(d=3)

10-1
SNR 0dB 10dB 20dB

T60=0.3 s

Fig. 2. Experimental results for 3 sources within 5 seconds.
to perform better than mapped MUSIC when these methods exploit
the same order statistics and there are sufficient observed data. The
reason for these results may be the statistical bias derived from theshortage of observed data, because cumulants tend to be affectedby statistical bias more than moments. Figure 5 shows the average
processing time under these conditions. As a result, our proposed
mapped MUSIC is sufficiently efficient in practical conditions witha realistic number of observed signals and has much less computational complexity than 2q-MUSIC.
5. CONCLUSION
In this paper, we proposed mapped MUSIC, a high-resolution DOA
estimator, which is applicable even in underdetermined conditions.
We described the mapped MUSIC algorithm and suitable mappingfor DOA estimation. We also discussed the property of the mapping function with the degree of d to analyze the 2d-th order crossmoments, and proposed efficient algorithms with which to calculate fourth and sixth order moments. Furthermore, we compared
the characteristics of the proposed method and the conventional
2q-MUSIC utilizing 2q-th order cumulants. We showed that the
proposed methods has much less computational complexity than
2q-MUSIC. We confirmed that DOA estimation with the proposed
mapped MUSIC technique is as accurate as that obtained withthe conventional 2q-MUSIC technique but with less computational
complexity.
6. REFERENCES

T60=0 sSNR 0dB 10dB 20dB

SNR 0dB 10dB 20dB

T60=0.12 s

T60=0.3 s

Fig. 4. Experimental results for 5 sources within 5 seconds.
103

Processing!!time![sec]RMSE [degree]

101

102
MUSIC
101

2q MUSIC!(q=2)
2q MUSIC!(q=3)

100Mapping!MUSIC!(d=2)
Mapping!MUSIC!(d=3)

10 1
10 2
3!sec

5 sec
DurationFig. 5. Processing time.
[3] P. Chevalier, A. Ferreol, and L. Albera, “High resolution direction finding from higher order statistics: The 2q-music algorithm,” IEEE Trans. Signal Process., vol. 54, no. 8, pp. 2986–
2997, 2006.
[4] H. Saruwatari, S. Kajita, K. Takeda, and F. Itakura, “Speechenhancement using nonlinear microphone array based on complementary beamforming,” IEICE Trans. Fundam., vol. E82-A,
no. 8, pp. 1501–1510, 1999.
[5] H. Kamiyanagida, H. Saruwatari, K. Takeda, F. Itakura, andK. Shikano, “Direction of arrival estimation using nonlinear
microphone array,” IEICE Trans. Fundam., vol. E84-A, no. 4,
pp. 999–1010, 2001.
[6] S. Miyabe, Biing-Hwang (Fred) Juang, Hiroshi Saruwatari, andKiyohiro Shikano, “Kernel-based nonlinear independent component analysis for underdetermined blind source separation,”in ICASSP, 2009, pp. 1641–1644.

[1] R. O. Schmidt, “Multiple emitter location and signal parameter
estimation,” IEEE Trans. Antennas Propag., vol. AP-34, no. 3,
pp. 276–280, 1986.

[7] J. A. Allen and D. A. Berkley, “Image method for efficiently
simulating small-room acoustics,” J. Acoust. Soc. Am., vol. 65,
no. 4, pp. 943–950, 1979.

[2] B. Porat and B. Friedlander, “Direction finding algorithms basedon higher order statistics,” IEEE Trans. Signal Process., vol. 39,
no. 9, pp. 2016–2024, 1991.

[8] I. A. McCowan and H. Bourland, “Microphone array post-filter
based on noise field coherence,” IEEE Trans. Speech Audio Process., vol. 11, no. 6, pp. 709–716, 2003.