DOA ESTIMATION FOR MULTIPLE SPARSE SOURCES
WITH NORMALIZED OBSERVATION VECTOR CLUSTERING
Shoko Araki†‡

Hiroshi Sawada†

Ryo Mukai†

Shoji Makino†‡† NTT Communication Science Laboratories, NTT Corporation
2-4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japan‡ Graduate School of Information Science and Technology, Hokkaido University
Kita 14, Nishi 9, Kita-ku, Sapporo-shi, Hokkaido 060-0814, Japan
Email: shoko@cslab.kecl.ntt.co.jp
ABSTRACTThis paper presents a new method for estimating the direction ofarrival (DOA) of source signals whose number N can exceed thenumber of sensors M . Subspace based methods, e.g., the MUSIC
algorithm, have been widely studied, however, they are only applicable when M > N . Another conventional independent component
analysis based method allows M ≥ N , however, it cannot be applied when M < N . By contrast, our new method can be applied
where the sources outnumber the sensors (i.e., an underdeterminedcase M < N ) by assuming source sparseness. Our method cancope with 2- or 3-dimensionally distributed sources with a 2- or 3dimensional sensor array. We obtained promising experimental results for 3 × 4, 3 × 5 and 4 × 5 (#sensors × #speech sources) in aroom (RT60 = 120 ms).
1. INTRODUCTION
Direction of arrival (DOA) estimation is a basic and important technique in array signal processing [1, 2]. The DOA estimation of speech,
on which we focus in this paper, has many applications such as teleconference systems and automatic conference minutes generators.
Such applications usually face situations where more sources thansensors are active simultaneously. We propose a new method for estimating DOA that can be applied even when the sources outnumber
the sensors.
The most widely used approaches for DOA estimation are subspace based methods, e.g., the MUSIC (MUltiple SIgnal Classification) algorithm [3], and its variants. Because these methods need anoise subspace, they require more sensors than sources M ≥ N + 1,
that is they can be applied only when M > N .
Recently, a DOA estimation method based on independent component analysis (ICA) has been proposed [4, 5]. It estimates DOAs
directly from the separation matrix estimated with ICA by utilizing
the fact that the separation matrix is related to the source mixing
process. Because this method is based on ICA, it allows cases where
M ≥ N . However, it still cannot be applied when M < N .
In order to cope with underdetermined cases where M < N , wepropose a new DOA estimation method by assuming source sparseness. If the signals are assumed to be sufficiently sparse in the timefrequency domain, we can assume that only one source is dominant
at a time-frequency point. Therefore, we can estimate information
about each source by collecting observation samples that appear tobelong to one of the sources. With our method, first we normalize
all the observations and cluster the normalized observation vectors.

1­4244­0469­X/06/$20.00 ©2006 IEEE

As each cluster corresponds to an individual source, we can formulate the DOAs by using the cluster centroids and given sensor locations. We have already proposed a sparseness based blind source
separation algorithm with normalized observation vector clustering[6]. In this paper, we show that we can also estimate the DOAs of
more sources than sensors with the same normalization and clustering method.
Some DOA estimation methods for M < N cases have alreadybeen proposed [7–9]. The method described in [7] is based on source
sparseness. It clusters the phase difference of only two sensor observations, then estimates the DOA. In addition, [8] estimated the DOA
with the cumulants of observations by assuming the non-gaussianityof sources, and utilizing a linear sensor array. These two methods
limit the DOA estimation ability on a 2-dimensional half-plane. The
authors of [9] utilized a triangular sensor array and clustered thephase differences of each sensor pair by assuming source sparseness. Their method expands the DOA estimation ability to an entire
2-dimensional plane. However, their approach still cannot handle 3dimensionally distributed sources. Moreover, their formulation assumed a regular-triangle sensor array. They have to re-formulate
their method to use an array with another arrangement.
On the other hand, our newly proposed formulation is moregeneral: the DOA is given by the cluster centroids and the sensor coordinates. Therefore, we can employ an arbitrary sensor arrangement including a 3-dimensional arrangement, and cope with
3-dimensionally distributed sources.
Experimental results show that the new method succeeded in estimating the DOAs for M × N of 3 × 4, 3 × 5 and 4 × 5. Neither
the MUSIC algorithm nor the ICA based method can be used in suchsituations. We also compare our method and MUSIC when M = 3
and N = 2. When sources were closely placed, the MUSIC algorithm failed to estimate their DOAs, whereas the proposed method
still succeeded.
2. PROBLEM DESCRIPTION
2.1. Observation model
Suppose that sources s1 , . . . , sN are convolutively mixed and observed at M sensors(1)
xj (t) = N
l hji (l) si (t − l), j = 1, . . . , M,
i=1
where hji (l) represents the impulse response from source i to sensor
j. In this paper, we especially consider a situation where the number of sources N can exceed the number of sensors M (M < N ).
Here we assume that the number of sources N is given or can beV ­ 33ICASSP 2006

(a)

qi
dj

qi
||qi||

OT

qi||qi||

qi||qi||

dJqi||qi||T(dj - dJ) qi||qi||

θjJi

(b) z

dj

x(f,τ)

qi

φi

y

θi

Sensorlocationsx

i=1 hi (f )si (f, τ ),≈hk (f )sk (f, τ ), k ∈ {1, · · · , N },(4)

where sk (f, τ ) is a dominant source at the time-frequency point(f, τ ). For instance this is true for speech signals in the time-frequency
domain [10, 11].
2.2. DOAs of source signals
Let us define the DOAs of sources more closely. Let dj and qi
be 3-dimensional vectors representing the locations of sensor j andsource i, respectively. We express the direction of source si as ||qqii ||with a unit-norm. Here, sensor locations dj are given.
In order to estimate the DOAs of sources, we approximate thefrequency response hji (f ) with a far-field model [Fig. 1(a)]:qihji (f ) ≈ exp [j2πf c−1 dTj],(5)||qi ||
where c is the propagation velocity of the signals. That is, we assume that the frequency response depends only on the path difference dTj ||qqii || between a source i to origin O and a source i to asensor j [Fig. 1(a)]. When we consider the two sensors j and J, we
obtain the following expressions:hji (f )qi≈ exp [j2πf c−1 (dj − dJ )T](6)
hJi (f )||qi ||=

exp [j2πf c−1 ||dj − dJ || cos θijJ ].

D

DOA estimation

qk||qk||

angle cos θijJ with respect to a sensor pair j-J, [see Fig. 1(a)]. We
term the former ||qqii || an absolute DOA, and the latter cos θijJ a relative DOA. The previous methods with two sensors [7] or a linear
sensor array [8] can estimate only the relative DOAs. On the other
hand, in this paper, our goal is to estimate the absolute DOAs ||qqii ||of sources from the all sensor observations x(f, τ ). Using the azimuth θi and elevation φi [see Fig. 1(b)], the absolute DOAs ||qqii ||can be written asqi(8)= [cos θi cos φi , sin θi cos φi , sin φi ]T .
||qi ||
3. PROPOSED METHOD(3)

where hji (f ) is the frequency response from source i to sensor j,si (f, τ ) is the STFT of a source signal si , and τ is a time index. We
call x = [x1 , . . . , xM ]T an observation vector and hi = [h1i , . . . ,hM i ]T is a vector of the frequency responses from source si to all
sensors.
In the time-frequency domain, the sparseness of a source signal,
which has been widely employed for solving the underdetermined
problem [6, 7, 9–11], becomes prominent, if the source is coloredand non-stationary such as speech. When the signals are sufficiently
sparse, we can assume that the sources rarely overlap at each timefrequency point, and (3) can be approximated asx(f, τ )ck

ways: the DOA ||qqii || with respect to a coordinate system, and theestimated. Our task is to estimate the directions of arrival (DOAs) ofN sources from the sensor observations without information aboutsources si or impulse responses hji . We will formulate the DOA
estimation problem in Sec. 2.2.
As with most DOA estimation techniques, this paper employs atime-frequency domain approach. Using a short-time Fourier transform (STFT), the convolutive mixtures (1) can be converted to instantaneous mixtures at each frequency f :
xj (f, τ ) ≈ N(2)
k=1 hji (f )si (f, τ ),N

Clustering

Fig. 2. Flow of proposed method

Fig. 1. (a) Far-field model, (b) definition of source location.

x(f, τ ) ≈x(f,τ)

dmax

O=[0,0,0]Tor in vector notation,

Normalize(7)

These two equations show that we can estimate the DOA in two

This section explains our proposed DOA estimation method, which
is applicable when M < N . When the sparseness assumption (4)holds, the observation vector x(f, τ ) is a multiplication by an unknown scalar sk of the frequency response hk (f ) about a source sk
which is dominant at (f, τ ). Moreover, remember that the frequency
response hk (f ) includes the DOA information ||qqkk || of source sk[see (5)]. Therefore, if we can collect the observation samples thatappear to belong to one of the sources, we can estimate its DOA as
an expectation value.
Based on this idea, we cluster the observation vector x(f, τ ) sothat each cluster corresponds to each source. Figure 2 shows theflow of our method. To make x(f, τ ) form clusters, first we neednormalizations of all the observations. Then we cluster the normalized observation vectors x̄(f, τ ) to gather the observation samples
that belong to one of the sources. Because each cluster corresponds
to an individual source, finally, we estimate the DOAs by using thecluster centroids ck . We explain each step in detail in the following
subsections.
3.1. Normalization
In this step, we normalize all observation vectors x(f, τ ) so that they
form clusters, each of which corresponds to an individual source.
First, we remove the inconstancy of the scalar sk (f, τ ) in (4).
This can be normalized by taking the ratio of two observation components: xj (f, τ )/xJ (f, τ ) ≈ hjk (f, τ )/hJk (f, τ ), which is modeled as (6). Moreover, if we can normalize the frequency dependence
observed in the phase component of (6), we can handle all the frequency components together and obtain enough data samples even if
we use short observations. The frequency dependence is normalized
by dividing the phase of xj (f, τ )/xJ (f, τ ) by f .
Including the above normalizations, our normalization is performed for the all components of x(f, τ ) by selecting one of thesensors J:
arg[xj (f, τ )/xJ (f, τ )]x̄j (f, τ ) ← |xj (f, τ )| exp j(9)4f c−1 dmax
where c is the propagation velocity and dmax is the maximum distance between the reference sensor J and a sensor ∀ j ∈ {1, . . . , M }.

V ­ 34

The rationale for the frequency normalization 4f c−1 dmax can befound in the Appendix. Here, we keep the amplitudes at all sensors,
although this is not needed for DOA estimation. This information isuseful when we perform source localization, i.e., estimate the positions qi of sources, using a near-field model.
We also employ unit-norm normalization to handle the vectors
on a unit-hypersphere,• Reverberation time: 120 ms• Room height: 240cm

for x̄(f, τ ) = [x̄1 (f, τ ), . . . , x̄M (f, τ )] .
If sources are sparse, each component of the normalized observation vector is expressed as


qk
1
π
x̄j (f, τ ) = √ exp j(dj − dJ )T
(11)
2dmax
||qk ||Mπ1||dj − dJ || cos θkjJ , (12)= √ exp j
2dmax
Mby using (4), (5), (9), and (10). We can see that the normalizedcomponents x̄j (f, τ ) keep the geometric information of sensors anda source sk , which is dominant at (f, τ ). Therefore, the normalized
vectors x̄(f, τ ) form clusters based on the source geometry in theclustering step.

ck ← E[x̄(f, τ )]k , ck ← ck /||ck ||,(13)

where E[·]k is a mean operator for the members of a cluster Ck . The
clustering criterion is to minimize the total sum J of the squared
distances between the cluster members and their centroid


2
J = M
(14)
x̄∈Ck ||x̄ − ck || .
k=1 Jk , Jk =This minimization can be performed efficiently with the k-meansclustering algorithm [12].
3.3. DOA estimation
Because each cluster corresponds to an individual source, the centroid ck represents the geometry of the source sk . Note that thecentroid gives an expectation value for the k-th source [see (13)].
From (11) and (12), the j-th component of the ck is expressed as{ck }j∝E[x̄j (f, τ )]k


π
1
q̃k
= √ exp j(dj − dJ )T
2dmax||q̃k ||Mπ1||dj − dJ || cos θ̃kjJ .
= √ exp j
2dmax
M

(15)
(16)

where ||q̃q̃kk || and θ̃kjJ are the estimated DOAs. When we use all thecomponents of the centroid ck , we can derive the following relationship from (15),
q̃k
2dmax
D
(17)
=
rk
||q̃k ||
π
where D = [d1 − dJ , · · · , dM − dJ ]T is a matrix of the sensor locations and r k = [arg[{ck }1 ], · · · , arg[{ck }M ]]. Using theMoore-Penrose pseudo-inverse, the absolute DOA is obtained asq̃k
2dmax +
(18)=
D rk .
||q̃k ||π
If rank(D) ≥ 3, we can estimate the 3-dimensional absolute DOA.
Note that if we use a unit sensor array system [d1 , d2 , d3 , d4 ] =α[(0, 0, 0)T , (1, 0, 0)T , (0, 1, 0)T , (0, 0, 1)T ] (α: a constant) and J =

282cm

x

S3

[Setup 2]
z
S4

50cm

S5

80cm

200cm

y

4cm

(10)

In this step, we find clusters C1 , . . . , CN formed by all normalized
vectors x̄(f, τ ). The centroid ck of a cluster Ck is calculated by

x

z
S1

T

3.2. Clustering[Setup 1]
y
4cm

y

S2

375 cmx̄(f, τ ) ← x̄(f, τ ) / ||x̄(f, τ )||880 cmx
4cm

Fig. 3. Experimental setups with non-linear sensor arrays.
Table 1. Experimental results for M = 3, N = 4 (Setup 1),
Sources1s2
s3s4
MUSIC (each)24◦ 117◦ 217◦ 311◦
Proposed (5 sec.) 25◦ 114◦ 214◦ 313◦
Proposed (1 sec.) 23◦ 112◦ 212◦ 318◦r
1, then arg[{ck }1 ] ≈ 0 and ||q̃q̃kk || = ||r k || , where r k = [arg[{ck }2 ,karg[{ck }3 ], arg[{ck }4 ]]T . That is, we do not need D + , and DOAq̃k
can be obtained simply from r k .
||q̃k ||
Using j-th component of the centroid ck , we can also estimate
the relative DOA θkjJ from (16) if needed:cos θ̃kjJ =2dmax arg [{ck }j ].
π ||dj − dJ ||(19)

4. EXPERIMENTAL RESULTS
We performed experiments in a reverberant condition. Observations
were made by following (1) with the impulse responses hji (l) measured in a room (Fig. 3) and 5-second English speech sources si (t)sampled at 8 kHz. The frame size L for STFT was 512, and theframe shift was 256 (= L/2).
We estimated the absolute DOAs ||qqkk || with (18), and evaluatedthem with the azimuth θk and the elevation φk [see (8)]. As the true
directions, we utilized the DOAs estimated by the MUSIC algorithm
when there was only one source signal. They will appear in thetables with the label “MUSIC (each)”. We investigated eight speakercombinations and averaged the results.
Table 1 shows the results for four sources with three sensors(M = 3, N = 4), that were arranged non-linearly (Fig. 3 [Setup1]). Here, all source heights were the same as the height of the sensor array. Because all elevations φk are zero, only the results of
θk are shown in the Table 1. We can see that the DOAs estimatedwith our proposed method were very close to the results of “MUSIC(each)”. Even when we used only one second data for the DOA estimation, we still obtained reasonable results as shown in Table 1. In a
more complicated situation where M = 3 and N = 5, the proposed
method estimated the DOAs very accurately as shown in Table 2.
We also applied our method to a 3-dimensional sensor arrangement (Fig. 3 [Setup 2]). In this case, the sources had different heights,
and therefore, we estimated both azimuths θk and elevations φk .
Table 3 shows results for five sources with four sensors (M = 4,
N = 5). As regards θk , although the estimation error was sometimes greater than the result for M = 3, N = 4 (Table 1), we
still obtained reasonable results for such a complicated case. TheV ­ 35

elevation values φk were also very close to the true values, “MUSIC (each)”. We can say that our proposed method can be applied
to such a 3-dimensional DOA estimation and that it gives us fairly
precise DOAs.
To show the effectiveness of our proposed method even for a situation where the MUSIC algorithm can be applied, we performedexperiments for two-source three-sensor (Setup 1) cases with bothmethods. Figure 4 shows the resolution of both methods, and Table 4shows the estimated DOAs θk . Figure 4 (a) and (c) are example MUSIC spectra at a frequency f = 1844 Hz, and (b) and (d) are DOA
histograms for members of each cluster (the DOA of each member
can be calculated with (18) using x̄(f, τ ) instead of the centroid ck ).
When two sources were placed with a wide spacing [Fig. 4 (a) and(b), and Table 4 “wide”], both MUSIC and the proposed method estimated the directions well enough. In contrast, when the two sourceswere close to each other [Fig. 4 (c) and (d), and Table 4 “close”],
MUSIC failed to estimate the two directions, whereas the proposed
method was still successful. We consider that our proposed methodwith sparseness assumption has a high resolution, although this resolution depends on the sparseness of the source signals and is affectedby the room reverberation condition.

−π/2 ≤ arg[x̄j (f, τ )] ≤ π/2

holds. Due to the frequency normalization with 4f c−1 dmax , thephase is distributed effectively from −π/2 to π/2. This property isalso important for the distance measure (14), since |x̄ − x̄ | increases
monotonically as | arg(x̄) − arg(x̄ )| increases.

0
0

120
240
DOA [deg.]

60

500

0
0

360

(c)

40

Histogram

500

(b)

120
240
DOA [deg.]

2000

360

(d)

1500
1000

20
0
0

120
240
DOA [deg.]

360

500
0
0

120
240
DOA [deg.]360

Table 4. Experimental results for M = 3, N = 2 (Setup 1),
wide
close
Sources1s2
s1s2
MUSIC (each)117◦ 217◦ 132◦ 154◦
Proposed (5 sec.) 114◦ 217◦ 128◦ 156◦
MUSIC (mix)125◦ 217◦ 126◦ 345◦
6. REFERENCES[1] S. U. Pillai, Array Signal Processing, Springer-Verlag, 1989.
[2] M. Brandstein and D. Ward, Eds., Microphone Arrays,
Springer, 2001.
[3] R. O. Schmidt, “Multiple emitter location and signal parameter
estimation,” IEEE Trans. Antennas and Propagation, vol. 34,
no. 3, pp. 276–280, Mar. 1986.
[4] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Frequencydomain blind source separation,” in Speech Enhancement,
J. Benesty, S. Makino, and J. Chen, Eds., pp. 299–327.
Springer, Mar. 2005.
[5] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Multiple
source localization using independent component analysis,” inProc. AP-S/URSI2005, 2005.
[6] S. Araki, H. Sawada, R. Mukai, and S. Makino, “A novel blindsource separation method with observation vector clustering,”
in Proc. IWAENC 2005, Sept. 2005, pp. 117–120.
[7] S. Rickard and F. Dietrich, “DOA estimation of many Wdisjoint orthogonal sources from two mixtures using DUET,”in Proc. SSAP2000, Aug. 2000, pp. 311–314.
[8] S. Shamsunder and G. B. Giannakis, “Modeling of nonGaussian array data using cumulants: DOA estimation of more
sources with less sensors,” Signal Processing, , no. 30, pp.
279–297, 1993.
[9] M. Matsuo, Y. Hioka, and N. Hamada, “Estimating DOA
of multiple speech signals by improved histogram mapping
method,” in Proc. IWAENC2005, Sept. 2005, pp. 129–132.

We proposed a new DOA estimation method for cases where M <N by assuming source sparseness. The method is based on the normalization and clustering of the observation vectors. We obtained
promising experimental results for underdetermined cases in a reverberant condition. We also confirmed that our proposed method
has higher resolution for estimating the directions of sources thanthe MUSIC algorithm.
We can also formulate the source localization, i.e., estimate thepositions qi of sources, using our proposed normalized observation
vector and a near-field model. We plan to investigate its performancein the future.

an inequality

1000

Fig. 4. Resolution of MUSIC [(a),(c)] and the proposed method[(b),(d)]. (a)(b): wide setting (θ1 =117◦ , θ2 =217◦ ), (c)(d): close setting (θ1 =132◦ , θ2 =154◦ ).

5. CONCLUSION

Appendix
This appendix explains the rationale for the denominator of thephase normalization in (9). The normalized observation vector component can be written as (11) or (12). From the fact that

 


 
jJ 
q
(dj − dJ )T ||qkk ||  = ||dj − dJ || cos θk  ≤ dmax ,

1000

(a)

MUSIC spectrum

Table 3. Experimental results for M = 4, N = 5 (Setup 2),
Sources1 s2 s3s4
s5
θ 31◦ 85◦ 133◦ 222◦ 302◦
MUSIC (each)φ −26◦ 6◦ 30◦ 39◦ −8◦θ 30◦ 79◦ 132◦ 221◦ 298◦
Proposed (5 sec.)
φ −22◦ 7◦ 28◦ 35◦ −9◦

1500

Histogram

MUSIC spectrum

Table 2. Experimental results for M = 3, N = 5 (Setup 1),
Source
s1 s2
s3
s4
s5
MUSIC (each) 24◦ 117◦ 176◦ 217◦ 311◦
Proposed (5 sec.) 23◦ 112◦ 175◦ 218◦ 314◦

[10] Ö. Yılmaz and S. Rickard, “Blind separation of speech mixtures via time-frequency masking,” IEEE Trans. SP, vol. 52,
no. 7, pp. 1830–1847, 2004.
[11] P. Bofill and M. Zibulevsky, “Blind separation of more sources
than mixtures using sparsity of their short-time Fourier transform,” in Proc. ICA2000, 2000, pp. 87–92.
[12] R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification,
Wiley Interscience, 2nd edition, 2000.

V ­ 36