A NEW SUBBAND INFORMATION FUSION METHOD FOR WIDEBAND DOA ESTIMATION

USING SPARSE SIGNAL REPRESENTATION

Ji-An Luo1,2 , Xiao-Ping Zhang1 , Zhi Wang2

1. Department of Electrical and Computer Engineering, Ryerson University

350 Victoria Street, Toronto, Canada M5B 2K3. Email: xzhang@ee.ryerson.ca
2. State Key Lab of Industrial Control Technology, Zhejiang University
Hangzhou, Zhejiang, China 310027. Email: {anjiluo, wangzhizju}@gmail.com

ABSTRACT signal subspaces by transforming the observation vectors as-

sociated with each bin into the focusing subspace and can deal
We present a new subband information fusion (SIF) method with coherent sources by averaging the subspace-aligned co-
for wideband direction-of-arrival (DOA) estimation using sin- variance matrices. Compared with ISSM, CSSM can enhance
gle sparse signal representation of multiple frequency-based DOA resolution and improve the accuracy of DOA estimates
measurement vectors. The problem of wideband DOA esti- at low SNR. However, CSSM requires an initial DOA esti-
mation using SIF method is to jointly utilize all the frequency mate and the precision of DOA pre-estimates greatly influ-
bin information to recover a single sparse indicative vector ences the accuracy of DOA estimation [4, 5].
(SIV). The SIF method belongs to the sparse signal represen-
tation domain and therefore it will suffer from two cases of Recently, a class of sparse signal representation (SSR)
ambiguity: algebraic aliasing and spatial aliasing. We show methods provide a new perspective for wideband DOA es-
that these two categories of ambiguity can be reduced by com- timation. The DOA estimation problem can be formulated
bining all the frequency components. The SIF algorithm is as recovering a spatial sparse signal vector or matrix by mini-
then proposed and the SIV is recovered iteratively. The nu- mizing the residual norm under sparsity constraint. One of the
merical simulations are performed to illustrate that the SIF most successive `1 -norm-based SSR algorithms for DOA esti-
method has superior performances. mation is `1 -SVD (Singular Value Decomposition) [6], which
reduces the computational complexity by SVD. Hyder and
Index Terms— Direction-of-arrival estimation, sparse Mahata [7] present a joint `2,0 -norm approximation (JLZA)
signal representation, subband information fusion, wideband method and extend it to wideband DOA estimation. The main
source, unconstrained optimization. limitation of the `1 -SVD and JLZA algorithms is that they
could not jointly use all the subband information to estimate
1. INTRODUCTION DOA and therefore lose their performance. More recently,
Liu et al. [8] present a wideband covariance matrix sparse
Many wideband direction-of-arrival (DOA) estimation meth- representation (W-CMSR) method for DOA estimation. The
ods have been proposed over last three decades due to their W-CMSR method uses time domain measurement informa-
various applications in radar, sonar, wireless communication tion and has its limitation for spatial nonambiguity because
and radio-astronomy etc. [1]. A classical wideband array the spatial aliasing is frequency-dependent. To deal with the
processing is to decompose the wideband signals into many spatial aliasing for SSR problem, Tang et al. [9] shows that
narrowband signals with a filter bank or the discrete Fourier such ambiguity can be removed by using multiple dictionar-
transform (DFT), and two categories, referred to as incoherent ies, each dictionary corresponding to a judiciously chosen fre-
signal subspace method (ISSM) [2] and coherent signal sub- quency. However, it still does not consider using all the fre-
space method (CSSM) [3], are utilized to realize wideband quency bin information to reduce the spatial ambiguity.
DOA estimation. The ISSM estimates the DOAs indepen- In this paper, we present a new subband information
dently and average them over all the bins. The performance fusion (SIF) method for wideband DOA estimation. This
of ISSM may deteriorate with low signal-to-noise ratio (SNR) method can jointly integrate all the frequency bins together
frequency bins and coherent sources. The CSSM align the to estimate a sparse indicative vector (SIV) that is used to
indicate the location of sources. By using the SIF method,
This work was supported in part by Canada NSERC Grant No. RG-
PIN239031, the National Natural Science Foundation of China (NSFC) un-
the spatial nonambiguity condition can be extended dramat-
der Grant No. 61273079, and the Strategic Priority Research Program of the ically compared with the classical beamforming technique.
Chinese Academy of Sciences under Grant No. XDA06020201. We then develop the SIF algorithm and compare the SIF al-

978-1-4799-0356-6/13/$31.00 ©2013 IEEE 4016 ICASSP 2013

gorithm with the W-CMSR algorithm. The simulation results where vk,n = [vk,1 (n), . . . , vk,L (n)]T is the sparse represen-
show that the proposed algorithm has better performance. tation of source vector. The nonzero entries of vk,n represent
true sources and zero otherwise. Assume that the nonzero in-
2. STATISTICAL MODEL dexes of vk,n are contained in the set I := {l1 , . . . , lQ }. Ob-
viously, vk,lq (n) = sk,q (n), where sk,q (n) denotes the k-th
Consider a uniform linear array (ULA) of P omnidirectional subband signal radiated from the q-th target at the n-th snap-
sensors is exposed to a wavefield generated by Q far-field shot.
wideband sources in presence of noise. For wideband pro-
cessing, a standard way is to split the time-samples at each 3. WIDEBAND DOA ESTIMATION USING A
sensor into N segments, where for each segment, K narrow- SUBBAND INFORMATION FUSION METHOD
band signals are obtained by a bank of narrowband filters or
the discrete Fourier transform (DFT) [5]. Assume that the 3.1. A New SIF Method
frequencies of all related sources overlap from ω0 − B2 rad/s
From the sparse representation model in (5), it is clear that
to ω0 + B2 rad/s, where ω0 is the center frequency and B is
the matrices V1 , . . . , VN share the identical sparse structure,
the bandwidth of some frequency band where the frequency
where Vn = [v1,n , . . . , vK,n ], n = 1, . . . , N . Let V =
bands of all sources intersect. The array output vector for
[V1 , . . . , VN ]. The nonzero rows of matrix V indicate the
fixed frequency ωk ∈ [ω0 − B2 , ω0 + B2 ], k = 1, . . . , K, is
source locations. The sparse indicative vector employed in
formed and can be expressed by
[10, 11] for narrowband DOA estimation can also be used to
xk,n = Ak (θ)sk,n + wk,n (1) represent the sparsity structure of V in wideband processing.
Thus, the estimating of entire matrix V can be avoided. To
where xk,n is the P × 1 array measurement vector, sk,n de- develop the SIF method, we introduce an over-complete dic-
notes the Q × 1 source signal vector at frequency ωk and wk,n tionary [10]
represents the P × 1 additive noise vector. Ak (θ) is the man-  
ifold matrix at frequency ωk Ψk,n = ak (θ1 )v̂k,1 (n) · · · ak (θL )v̂k,L (n) , (6)

Ak (θ) = [ak (θ1 ), . . . , ak (θQ )] (2) where

−1
where the steering vector ak (θq ) can be written as v̂k,l (n) = aH H
k (θl )ak (θl ) aH
k (θl )xk,n . (7)
iT
Let us use D ∈ CL×L denote a new dictionary that DH D =
h d d
ak (θq ) = 1, e−jωk c sin θq , . . . , e−jωk (P −1) c sin θq (3)
Φ, where
K X N
where θq , q = 1, . . . , Q, is the DOA of the q-th source, d is
X
Φ= ΨHk,n Ψk,n . (8)
the distance between the two adjacent sensors, and c is the k=1 n=1
speed of the signal propagation.
The DOA estimation appears in (1) is a nonlinear pa- One realization of D is given by the eigenvalue decomposition
rameter estimation problem, where the DOA parameter θ = (EVD) of the Hermitian matrix Φ. Assume that Φ can be
[θ1 , . . . , θQ ]T need to be estimated. The sparse representation written as Φ = UΛUH , where U is a unitary matrix whose
model transforms a parameter estimation problem into sparse columns are composed by L orthonormal eigenvectors and Λ
L
spectrum estimation. We denote Θ = {θl }l=1 as the set of a is a diagonal matrix of eigenvalues. Thus, the dictionary D
1

sampling grid of all possible source locations, L  Q. We has the expression D = UΛ 2 . Let M denote the rank of Φ.
assume that the grid is fine enough that Θ can represent the We have rank(D) = rank(Φ) = M .
true source locations exactly or closely. Then, the array out- We define a new SIV g which has the same sparsity
put model (1) is changed into structure with V. g can be used to indicate the locations of
sources. Recovering g needs to jointly minimize a residual
2
L
X term ky − Dgk and a `p -norm (0 ≤ p ≤ 1) penalty, which
xk,n = ak (θl )vk,l (n) + wk,n (4) leads to the SIF method:
l=1
2
minf (g) : f (g) = ky − Dgk + λkgkp , (9)
where vk,l (n) and ak (θl ) are the visual source and steer- g

ing vector corresponding to the l-th grid. We introduce an

where λ is a regularization parameter, y is a new vector de-
overcomplete basis matrix Ak = [ak (θ1 ), . . . , ak (θL )]. The
fined by DH y = h,
sparse representation model in (4) can be represented by the
compact model K X
X N
h= ΨH
k,n xk,n . (10)
xk,n = Ak vk,n + wk,n (5) k=1 n=1

4017
If the number of sources is fixed during N snapshots, g is Q- use the unique solution condition (see [13], theorem 2.4) for
sparse and kgk0 = Q. The indexes of nonzero elements in g the SMV problem. At the same time, the unique solution for
are contained in the set I which indicates the source locations. y = Dg is under the MMV problem, since it is derived from
Intuitively, g captures the sparsity structure of y in the basis the MMV problem.
of D. Therefore, g can be recovered as sparse as possible We then investigate the spatial aliasing problem. Assume
provided that the residual term ky − Dgk2 is minimized. The that two angles θl1 and θl2 satisfy
problem of DOA estimation can then be solved by recovering
g instead of estimating V. ωk d sin θl1 /c = ωk d sin θl2 /c + 2πI, (12)
Note that the derivation of (9) has an equivalent expres- where I is an arbitrary integer. Let us use θl1 and θl2 to gener-
sion which is given by ate two corresponding columns in dictionary D. If these two
columns are identical, the spatial aliasing occurs. Unlike the
arg minf (g) narrowband DOA estimation problem, the spatial nonambi-
g
K X
X N guity condition for wideband DOA estimation can be more
≡ arg min kxk,n − Ψk,n gk2 + λkgkp . (11) relaxed. It is clear that the spatial aliasing is dependent on the
g
k=1 n=1 frequency. For the SIF method, all the frequency bin informa-
tion can be combined to reduce the spatial ambiguity, which
The above expression (11) shows that all the frequency bins leads to the following theorem.
and snapshots are integrated together to estimate g. We then
find that the over-complete dictionary Ψk,n can only be used Theorem 2. Assume that the number of sources satisfies the
for specific frequency and snapshot in which xk,n can be ex- unique solution condition proposed in theorem 1. Let ∆ω =
pressed as a sparse signal. Yet D includes all the frequency ωk − ωk−1 . With a ULA interspaced by a unit distance d, the
πc
and snapshot measurements information and therefore (9) is spatial nonambiguity is guaranteed if d < ∆ω , where c is the
actually an information fusion formulation for the wideband propagation speed of the source signal.
DOA estimation problem. Proof : See [12]. The proof of theorem 2 is similar to the
proof of theorem 1 in [9]. We first assume that θl1 is one of
3.2. Nonambiguity Guarantee for the SIF Method the target angles and θl2 is the spatial aliasing angle, which
lead to the identical column as θl1 in Hermitian matrix Φ. We
As discussed in [9], the sparse signal representation (SSR) then prove that under the condition d < ∆ω πc
, θl1 = θl2 holds
approaches are generally subject to two kinds of ambiguity. and the spatial nonambiguity is guaranteed.
One is referred as algebraic aliasing, which arises from the
over-completeness of the dictionary. The other is called spa-
3.3. SIF Algorithm
tial aliasing, just like the classical beamforming. If the array
spacing d is larger than half the apparent wavelength, it will Consider the following unconstrained problem:
be possible to find at least two angles θl and θl0 . The corre-
K X
N
sponding columns in Ak will be identical. X
min kxk,n − Ψk,n gk2 + λkgkp . (13)
For the algebraic aliasing-free problem, we investigate the g
k=1 n=1
sufficient condition for the existence of a unique solution to
the noise-free equations y = Dg. The following theorem The problem (13) can be solved by generating a sequence
shows that the algebraic nonambiguity is guaranteed if the of iterates {g0 , g1 , . . .}. For each iteration, the optimization
number of sources satisfies certain condition. problem of (13) is given by

Theorem 1. Consider the equations y = Dg, where D is a gt+1 = 2Π(gt )H−1 (gt )h, (14)
L × L dictionary, y is a L × 1 vector and g is a sparse in- where
dicative vector whose nonzero elements indicate the locations H(g) = 2ΦΠ(g) + pλI (15)
of sources. D and y are defined by DH D = Φ, DH y = h
(g12 1−p/2
 
respectively, where Φ and h is given in (8) and (10). As- + ) 0
sume that no spatial nonambiguity exists. A unique solution of Π(g) = 
 ..  , (16)

.
sparse indicative vector g is guaranteed if Q ≤ min(d(M + 0 2
(gL + )1−p/2
1)/2e − 1, d(P + N K)/2e − 1) for N K < P and Q ≤
min(d(M + 1)/2e − 1, P − 1) for N K ≥ P , where d·e de- I is a L×L identity matrix,  is a very small smooth parameter
notes the ceiling operation. and the expressions of Φ and h are given in (8) and (10). The
computational complexity of computing (14) can be reduced
Proof : See [12]. In Subsection 3.1, we convert the mul- by using the conjugate gradient (CG) technique [14, 15]. Let
tiple measurement vectors (MMV) problem to a single mea- H−1 (g)h be replaced by cg(H(g), h), where cg(H(g), h) is a
surement vector (SMV) problem. The proof of theorem 1 can CG solution of linear equation H(g)b = h for b = H−1 (g)h.

4018
 0
4. EXPERIMENTAL RESULTS
W−CMSR
−20 SIF

In this section, we illustrate the performance of SIF algorithm −40

via various numerical simulations. We consider the same nu-
−60
merical examples given in [8] for ease of comparison. As-

Power (dB)
sume that two BPSK signals with central frequency of 70 −80

MHz and bandwidth of 40% impinge on a ULA with 7 sen-

−100
sors. Fig. 1 and Fig. 2 show the DOA estimation results of two
sources from the directions of −10◦ and 10◦ obtained by the −120

W-CMSR and SIF respectively. We take K = 256, N = 1 −140

and SNR = 0 dB. The parameters p, λ and  are set to 0.1,
λ = 0.2 × k2hk∞ and  = 10−18 for the SIF algorithm.  −160

can be arbitrary small. In Fig. 1, the interspace d of the ULA −180

satisfies d = πc/(ω0 + B2 ). We then extend d to 100 times −80 −60 −40 −20 0
DOA (degree)
20 40 60 80

in Fig. 2. It can be seen from Fig. 1 that both the W-CMSR

and SIF algorithms perform well. However, the W-CMSR al- Fig. 2. Spatial spectrum of two sources from −10◦ and 10◦
gorithm fails in Fig. 2 while the SIF algorithm still has good obtained by the W-CMSR and SIF algorithms, P : 7, SNR: 0
results. Fig. 2 shows that the spatial ambiguity is reduced by dB, grid resolution: 1◦ , d: 194.8.
the SIF algorithm.
0
SIF:d=1.948
0 −20 SIF:d=194.8
W−CMSR
−20 SIF −40

−40 −60
Power (dB)

−60 −80
Power (dB)

−80 −100

−100 −120

−120 −140

−140 −160

−160 −180
−80 −60 −40 −20 0 20 40 60 80
DOA (degree)
−180
−80 −60 −40 −20 0 20 40 60 80
DOA (degree)
Fig. 3. Spatial spectrum of two sources from 6◦ and 10◦ ob-
tained by the SIF algorithm, P : 7, SNR: 0 dB, resolution 1◦ .
Fig. 1. Spatial spectrum of two sources from −10◦ and 10◦
obtained by the W-CMSR and SIF algorithms, P : 7, SNR: 0
dB, grid resolution: 1◦ , d: 1.948. 5. CONCLUSION

In this paper, we present a new subband information fusion

The interval of two adjacent frequencies is ∆ω = 2π in (SIF) algorithm to solve the wideband DOA estimation prob-
the simulations. According to the theorem 2, the spatial non- lem. The SIF algorithm utilizes all the frequency bin in-
ambiguity is guaranteed if d < 2c , c = 3 × 108 is the propa- formation to recover a sparse indicative vector (SIV) itera-
gation speed of the signal. As d increases, the SIF algorithm tively. We show that the algebraic ambiguity resulting from
captures the unique solution of (13). Moreover, the increas- the over-complete dictionary can be alleviated by multiple
ing of d makes the columns of dictionary D more and more measurement vectors. Since the spatial aliasing is frequency-
uncorrelated. In Fig. 3, the dash line plots the“mismatch” of dependent, the spatial ambiguity due to spatial aliasing can
the DOA estimate when two sources are close at 4◦ and 10◦ . be reduced by using all the frequency information. By using
The DOA estimation result for the dash line in Fig. 3 is 2◦ and the CG method, the computational complexity of the SIF al-
12◦ when d = 1.948. If the interspace d is increased to 100 gorithm has order O(L2 ). The performance of the proposed
times, the DOA estimation result exactly matches the source algorithm is demonstrated by numerical simulations.
directions, see Fig. 3 the solid line.

4019
 6. REFERENCES [12] J. A. Luo, X.-P. Zhang, and Z. Wang, “A new subband
information fusion method for wideband DOA estima-
[1] H. Krim and M. Viberg, “Two decades of array signal tion using sparse signal representation,” IEEE Transac-
processing research: the parametric approach,” IEEE tions on Signal Processing, prepared to be submitted.
Signal Processing Magazine, vol. 13, no. 4, pp. 67–94,
1996. [13] J. Chen and X. Huo, “Theoretical results on sparse rep-
resentations of multiple-measurement vectors,” IEEE
[2] G. Su and M. Morf, “The signal subspace approach for Transactions on Signal Processing, vol. 54, no. 12, pp.
multiple wide-band emitter location,” IEEE Transac- 4634–4643, Dec. 2006.
tions on Acoustics, Speech and Signal Processing, vol.
31, no. 6, pp. 1502–1522, 1983. [14] M.R. Hestenes and E. Stiefel, “Methods of conjugate
gradients for solving linear systems,” Journal of Re-
[3] H. Wang and M. Kaveh, “Coherent signal-subspace pro- search of the National Bureau of Standards, vol. 49, pp.
cessing for the detection and estimation of angles of ar- 409–436, 1952.
rival of multiple wide-band sources,” IEEE Transac-
[15] Z. He, A. Cichocki, R. Zdunek, and S. Xie, “Improved
tions on Acoustics, Speech and Signal Processing, vol.
FOCUSS method with conjugate gradient iterations,”
33, no. 4, pp. 823–831, 1985.
IEEE Transactions on Signal Processing, vol. 57, no.
[4] E. D. Di Claudio and R. Parisi, “WAVES: weighted av- 1, pp. 399–404, 2009.
erage of signal subspaces for robust wideband direction
finding,” IEEE Transactions on Signal Processing, vol.
49, no. 10, pp. 2179–2191, 2001.

[5] Y. S. Yoon, L. M. Kaplan, and J. H. McClellan, “TOPS:

new DOA estimator for wideband signals,” IEEE Trans-
actions on Signal Processing, vol. 54, no. 6, pp. 1977–
1989, 2006.

[6] D. Malioutov, M. Cetin, and A.S. Willsky, “A sparse

signal reconstruction perspective for source localization
with sensor arrays,” IEEE Transactions on Signal Pro-
cessing, vol. 53, no. 8, pp. 3010–3022, Aug. 2005.

[7] M.M. Hyder and K. Mahata, “Direction-of-arrival esti-

mation using a mixed `2,0 norm approximation,” IEEE
Transactions on Signal Processing, vol. 58, no. 9, pp.
4646–4655, Sept. 2010.

[8] Z. M. Liu, Z. T. Huang, and Y. Y. Zhou, “Direction-

of-arrival estimation of wideband signals via covariance
matrix sparse representation,” IEEE Transactions on
Signal Processing, vol. 59, no. 9, pp. 4256–4270, 2011.

[9] Z. Tang, G. Blacquiere, and G. Leus, “Aliasing-free

wideband beamforming using sparse signal representa-
tion,” IEEE Transactions on Signal Processing, vol. 59,
no. 7, pp. 3464–3469, 2011.

[10] J. A. Luo, X.-P. Zhang, and Z. Wang, “Direction-of-

arrival estimation using sparse variable projection op-
timization,” in Proc. Int. Symposium on Circuits and
Systems (ISCAS), 2012, pp. 3106–3109.

[11] J. A. Luo, X.-P. Zhang, Z. Wang, and M. Bao, “A new

direction-of-arrival estimation method based on the reg-
ularized sparse variable projection optimization,” IEEE
Transactions on Signal Processing, submitted.

4020