Direction-of-Arrival Estimation of Wideband Signals via ... - IEEE Xplore

64 downloads 0 Views 1MB Size Report
Abstract—This paper focuses on direction-of-arrival (DOA) estimation of wideband signals, and a method named wideband covariance matrix sparse ...
4256

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

Direction-of-Arrival Estimation of Wideband Signals via Covariance Matrix Sparse Representation Zhang-Meng Liu, Student Member, IEEE, Zhi-Tao Huang, and Yi-Yu Zhou

Abstract—This paper focuses on direction-of-arrival (DOA) estimation of wideband signals, and a method named wideband covariance matrix sparse representation (W-CMSR) is proposed. In W-CMSR, the lower left triangular elements of the covariance matrix are aligned to form a new measurement vector, and DOA estimation is then realized by representing this vector on an over-complete dictionary under the constraint of sparsity. The a priori information of the incident signal number is not needed in W-CMSR, and no spectral decomposition or focusing is introduced. Simulation results demonstrate the satisfying performance of W-CMSR in wideband DOA estimation in various settings. Moreover, theoretical analysis and numerical examples show how many simultaneous signals can be separated by W-CMSR on typical array geometries, and that the half-wavelength spacing restriction in avoiding ambiguity can be relaxed from the highest to the lowest frequency of the incident wideband signals. Index Terms—Convex optimization, direction-of-arrival (DOA) estimation, sparse representation, wideband signal.

I. INTRODUCTION

W

IDEBAND signals are widely used in various radar and sonar systems, and many methods have been proposed to estimate their directions-of-arrival (DOA), such as ISSM [1] and CSSM [2]. Most of those methods decompose the incident wideband signals into narrowband components, and then realize wideband DOA estimation with incoherent or coherent techniques. The coherent spectral focusing-based methods were demonstrated to surpass their incoherent counterparts [2], thus have attracted more research interest. However, there are two significant disadvantages within this kind of methods. First, DOA pre-estimates of the incident signals are required for spectral focusing, and the precision of those pre-estimates largely influences the performance of DOA estimation [3]–[5]. Second, they need the a priori information of the incident signal number,

Manuscript received November 30, 2010; revised March 23, 2011; accepted May 23, 2011. Date of publication June 09, 2011; date of current version August 10, 2011. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Jean Pierre Delmas. This work was supported in part by the National Natural Science Foundation under Grant 61072120 and 60802053, the Program for New Century Excellent Talents in University (NCET), Hunan Provincial Innovation Foundation for Postgraduate and the Innovation Foundation for Outstanding Postgraduates in the National University of Defense Technology under Grant B090402. The authors are with the College of Electronic Science and Engineering, National University of Defense Technology, Changsha, Hunan, China. (e-mail: [email protected]; [email protected]). This paper has supplementary downloadable multimedia material available at http://ieeexplore.ieee.org provided by the authors. This includes matlab codes. rar, which shows the source codes of W-CMSR. This material is 6.42 KB in size. Digital Object Identifier 10.1109/TSP.2011.2159214

which may not be available, especially in the noncooperative scenarios. Therefore, substitutive methods are required to explore better solutions for the problem of wideband DOA estimation. The technique of sparse representation provides a new perspective for DOA estimation. Methods of this category recover the spatial distribution of the incident signals by directly representing the array output on an over-complete dictionary under sparsity constraint, and the prior information of signal number is no longer a necessity. By now, several literatures have directly addressed sparsity-based DOA estimation [6]–[8]. The first one, as far as we know, is the global matched filter (GMF) method proposed by Fuchs [6]. GMF bases on uniform circular arrays, and exploits the beamformer samples to realize DOA estimation. Then Malioutov et al. proposed the method of L1-SVD [7], which reduces the dimension of the observations via singular decomposition. L1-norm-based optimization model achieved great success in those two methods. But to be more rigorous, the sparsity constraint on the spatial distribution of the incident signals should be modeled with L0-norm, so that the model parsimony is directly concerned by the sparsity constraint. However, the L0-norm-based optimization problem is not globally convergent, and it is NP-hard to solve it directly [9], so convex approximation is introduced and L1-norm is used in the above two methods to facilitate the optimization process [10], [11]. More recently, Hyder and Mahata introduced their L0-norm-based joint sparse approximation technique [12] into DOA estimation. In their work, the L0-norm constraint is approached by a family of convex functions, and a method called JLZA-DOA is proposed [8]. Although satisfying performance was reported for JLZA-DOA, global convergence is not guaranteed. L1-SVD and JLZA-DOA have also been extended to wideband signals. However, wideband L1-SVD and JLZA-DOA decompose the array output into narrowband components and obtain spectra in discrete frequency bins, so further efforts are required to derive integrative wideband DOA estimates from them. In this paper, we address the wideband DOA estimation problem by representing the array output covariance matrix under the constraint of sparsity. A method called wideband covariance matrix sparse representation (W-CMSR) will be proposed. In W-CMSR, we first align the lower left triangular elements of the array output covariance matrix to form a new measurement vector, and then reconstruct this vector on an over-complete dictionary to realize wideband DOA estimation.1 Due to the representation model that we use, the a priori information of source number is not required in W-CMSR. Meanwhile, W-CMSR processes the incident wideband signals holistically 1We have accomplished covariance matrix-based narrowband DOA estimation in a previous work [13], however, it was effortful to introduce this idea to wideband signals, since the extension is not so straightforward as it appears.

1053-587X/$26.00 © 2011 IEEE

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

and do not decompose them in the frequency domain, thus it does not get into the trouble of integrating narrowband signal components to obtain the final wideband DOA estimates. Further analysis shows that W-CMSR owns an intrinsic superiority in making good use of the array geometry to enhance its ability of separating multiple simultaneous signals, which may achieve the inspiring watershed of more signals than sensors in well-designed array geometries, such as the minimum redundancy array (MRA) [14]. Moreover, since the signal information within the signal band is exploited holistically, the half-wavelength restriction can be relaxed for ambiguity avoidance. This paper contributes to the area of wideband DOA estimation in the following five aspects. 1) W-CMSR relies less on the a priori information of the incident signal number than the ordinary subspace-based methods. 2) No spectral decomposition or focusing is introduced, thus W-CMSR is immune to imperfect DOA pre-estimates, and will not run into the tremendous efforts of focusing matrix selection or nonidentical narrowband DOA estimate fusion. 3) The a priori information of signal spectrum can be exploited in W-CMSR to improve the performance of DOA estimation. 4) Well-designed array geometries help W-CMSR to obtain enhanced ability in separating even more simultaneous signals than sensors. 5) The half-wavelength spacing restriction in avoiding ambiguity is relaxed from the highest to the lowest signal frequency, just like what JLZA-DOA has achieved. The rest of the paper mainly consists of six sections. Section II reviews the model of wideband DOA estimation. Section III presents the method of W-CMSR and its simplified variation when repeated off-diagonal elements exist in the perturbation-free covariance matrix. Section IV discusses how many signals can be separated in special array geometries and how much the half-wavelength restriction can be relaxed. Section V introduces a way for the fitting error threshold computation during L1-norm-based sparsity-inducing optimization in W-CMSR. Section VI carries out simulations to demonstrate the performance and computational efficiency of W-CMSR. Section VII concludes the whole paper. II. PROBLEM FORMULATION Suppose that wideband signals impinge onto an -element array from the directions of , respec-

4257

snapshots are collected, the snapshot at time is

tively, and given by

(1) is the waveform of the signal, is the propwhere agation delay of the signal from the reference point to the sensor, is the additive noise at the sensor, the time instants have been normalized with respect to the sampling interval to simplify notation. In this paper, we further assume that the incident signals and the additive noise are all Gaussian distributed and mutually independent. The waveforms of the independent signals satisfy and , with being the power of the signal, and being its unified correlation function defined as that satisfies (we unify it to make this correlation funcdenotes the unit impulse tion immune of the signal power), function, and is the conjugate operator. The unified correlation functions of the incident signals are assumed to be identical, i.e., , . This assumption tallies with many practical applications, such as digital communications where the signals are of similar modulation [15]. Based on this assumption, we simplify the unified correlation function of all the incident signals to . Further denote the corresponding power spectrum density function by , the effective bandwidth by . The number of snapshots is assumed to satisfy , with being the sampling interval. The propagation delay of the incident signals along the array in (1) can not be simplified to phase-delays due to the significant bandwidth, so spectral decomposition is always introduced before the narrowband DOA estimation methods can be applied. But during the procedure of model transformation, some of the signal energy may be lost [3], [5], and the a priori information of signal number is required to complete the narrowband DOA estimation process [16]. III. WIDEBAND DOA ESTIMATION VIA COVARIANCE MATRIX SPARSE REPRESENTATION The wideband DOA estimation method to be presented in this section bases on the correlation functions of the wideband signals, which can be derived from the array output covariance matrix. The perturbation-free covariance matrix of the wideband array output is given by (2), shown at the bottom of the page,

(2) .. .

.. .

..

.

.. .

4258

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

where is the variance of the additive noise. The directions of the incident signals can be derived from the correlation family for each if they are sepof arable, since there is a one-to-one mapping between the signal direction and the correlation family if no ambiguity occurs. But correlations are mixed in each covariance element, as those and the diagonal covariance elements are contaminated by the unknown noise variance, it is not an easy task to estimate the signal directions from directly. In this paper, we propose to use the correlation function to bridge the gap between them. In this section, we first present the wideband DOA estimation method, named W-CMSR, based on general array geometries, and then seek to save its computational load when repeated offdiagonal elements exist in . A. General W-CMSR The covariance matrix of the array output is conjugate symmetric, thus its upper right triangular elements can be represented by the lower left triangular ones. And as the main diagonal elements are contaminated by the unknown noise variance, we slide over them and align the lower left triangular elements column-by-column to obtain the following measurement vector de[see (3), shown at the bottom of the page], where notes the element of . Obviously, this vector can be decomposed into components as

sparse vector with nonzero values indexed corresponding to the signal direction locations in . Based on the above model, the following constrained sparsity-enforcing objective function can be introduced for wideband DOA estimation,

where

is

a

(8) denotes the L0-norm and indicates the spatial where distribution of the incident signals. The formulation of the over-complete dictionary plays an important role during the DOA estimation process based on (8). According to (6), the dictionary atoms depend on the unified correlation function of the incident signals. Two kinds of techniques can be used for the estimation of this function. One exploits the a priori information of the signal spectrum, while the other estimates it from the array output. In some scenarios, the correlation function can be computed from the prior information of signal modulation. For example, the amplitude of the unified correlation functions of phase-shiftkeying (PSK) signals own a triangular shape approximately, and (0,1) [15], thus if with its three vertexes located at the code rate is given or estimated with other methods, the unified correlation function can be computed straightforwardly as follows: (9)

(4) where is the carrier frequency. If no such prior information is available, we can first estimate the signal power spectrum from the array output, and then compute the unified correlation function as follows [17]:

with

(5)

(10)

Those signal components depend on the unified correlation function of the incident signals and their propagation delays along the array, which are directly related to the signal directions. Thus, the signal directions can be estimated from if the signal components can be separated. Denote the propagation delay of a signal at direction from the reference point to the sensor by , then the corresponding signal component with unit power is given by

where and the integral scope is set with respect to the signal bandwidth. After the formulation of , we introduce convex approximation to facilitate the solving of the NP-hard L0-norm-based optimization problem in (8). By further taking into consideration the perturbation contaminated in the practical measurement vector , we transform (8) to the following fitting error-constrained L1-norm-based optimization problem [10], [18], [19]: (11)

(6) Sample the potential spatial scope of the incident signals to get a direction set , e.g., if one discretizes the scope with interval of , then , thus can be rewritten in an over-complete form as (7)

where is the estimate of and is the fitting error threshold. The feasibility and efficiency of (11) for DOA estimation depends heavily on a reasonable selection of . In this paper, we set this threshold according to the perturbation level of . In practice, the array output covariance matrix is estimated from the snapshots, (12)

(3)

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

and the corresponding measurement vector is given by

4259

where (13)

where denotes the ance of the perturbation within

element of . The variis (14)

and is the propagation delay of a signal from direction between adjacent sensors. The perturbation variance within can be proved to satisfy

With straightforward derivation given in Appendix A, we conclude that

(19)

(15) where Thus we set the fitting error threshold in (11) in a heuristic way as follows to guarantee certain detection probability [20],

(16) where

is a weighting factor.

B. Simplified W-CMSR in Special Arrays When more than one pair of the array sensors are equallyspaced, repeated elements exist in . Those elements increase the problem dimension without providing any innovative information. Therefore, we propose combining them to simplify W-CMSR and save computation load. To facilitate notation, we take the uniform linear array (ULA) for example, but the analysis in this subsection can be extended to other array geometries alike with appropriate modifications. The perturbation-free covariance matrix of the ULA has a conjugate symmetric Toeplitz structure, and thus can be represented by its elements along the first column. Denote the propagation delay of the signal between adjacent sensors by , . After averaging the lower left covarithen ance elements diagonally, we obtain the following measurement vector: (17) . Then the spatial where distribution of the incident signals can be estimated by solving the following constrained L1-norm-based optimization problem [10], [18], [19]: (18)

where is defined the same as that in (15). Proof: See Appendix B. Thus, the fitting error threshold in (18) can be set as shown in (20) at the bottom of the page [20], where is a weighting factor. We delay the computation of the perturbation variances to Section V, and the selection of the weighting factors and to Section VI, so that we can concentrate on the main ideas of W-CMSR. In order to distinguish the two implementations of W-CMSR in the above two subsections, we denote the one in Part A by W-CMSR(GLA), with the identifier (GLA) (GLA: general linear array) introduced to indicate that this implementation bases on general array geometries, and possibly repeated elements in is not considered. We then denote the implementation in Part B by W-CMSR(ULA), with the identifier (ULA) used to indicate the incorporation of the repeated elements in based on special array geometries, such as ULA, but it should be clarified that this implementation also adapts to other arrays with pairs of equally-spaced sensors. After selecting the fitting error thresholds according to (16) and (20), we need to solve the constrained L1-norm optimization problems, i.e., (11) and (18), to finally obtain DOA estimates. In order to focus on the problem of wideband DOA estimation, we do not go any further into the detailed algorithms for solving such problems in this paper, and refer the interested readers to related literatures, e.g., [10], [18], [19], and the references therein, for more thorough discussions. IV. COMPLEMENTARY ISSUES TO W-CMSR In this section, we discuss some complementary issues to W-CMSR, including the separable signal number and halfwavelength relaxation. The discussions in this section base on and , but the perturbation-free assumption, i.e., the conclusions also provide constructive guidelines for the behavior of W-CMSR in practical applications.

(20)

4260

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

A. Separable Signal Number

Thus

Based on Corollary 1.2 of [21], we have demonstrated how many signals can be separated by narrowband CMSR via numerical examples in [13]. Simulation results show that by performing the sparse representation process on the covariance vector in the narrowband settings, one can make good use of the array geometry to separate more signals than sensors in well-designed arrays, such as MRA. We extend this conclusion to wideband signals here and present the following theorem on the separable signal number for W-CMSR. Theorem 1: With linear arrays interspaced by multiples of a unit distance that satisfy the nonambiguity constraint within the signal band, and the distances between pairs of sensors form , with being the array aperture, the set then W-CMSR can separate up to simultaneous wideband sig, and nals. Particularly, this number in -element ULA is that in MRA may exceed the sensor number. Proof: Select the pairs of sensors spaced by , and align the correlations of their outputs to form the measurement vector , then is a sufficient statistic signals impinge simultaneously of in (3). Assume that onto the array, denote the propagation delay of the signal along the distance of by , then can be decomposed into signal components as follows:

(21) where As

. , (21) can be rewritten as

(22) with (23), shown at the bottom of the page, where , , and Suppose that (21) also holds for another signals, with their powers satisfying . Deby note their propagation delays along a distance of , then

(24)

(25) Denote the real and imaginary parts of by and , those of by and , then , , and . As , and are all real-valued, thus splitting (25) into real and imaginary parts and aligning them column-wise yields

(26) Define

then is a continuous and bounded function with respect to . And for most commonly used signals, such as the PSK ones, and is continuous for . Thus, is also a continuous function. Then it can be derived from within satisfying (26) that, there is a certain frequency . As , one can further conclude that . Synthesizing the real and imaginary parts of yields (27) to spatial moments of inciwhich indicates that the dent signals at frequency equal those of another signals . According to Corollary 1.2 correspondingly, and , thus proves in [21], this proposition never holds when Theorem 1 by contradiction. B. Half-Wavelength Relaxation for Nonambiguity Guarantee In the method of JLZA-DOA, the signal components of different frequencies within the signal band are recovered jointly, thus the half-wavelength nonambiguity constraint is relaxed to the lowest signal frequency [8]. The temporal correlation-based method of W-CMSR also has an inherent property of using the wideband signal spectrum holistically, so it is also expected to relax the half-wavelength restriction in ambiguity avoidance. This property of W-CMSR is given in Theorem 2.

(23)

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

Theorem 2: If half-wavelength nonambiguity is guaranteed at the lowest signal frequency, it is holistically guaranteed in W-CMSR for wideband DOA estimation. Proof: Suppose that half-wavelength nonambiguity in DOA estimation at the lowest signal frequency is satisfied for the array, we divide the signal bandwidth into the higher and and , according to the criterion lower parts, denoted by of potential ambiguity, then . Suppose that another signals are ambiguous to the true signals within . According to (24), we have

4261

(32)

(28)

relies on the unified correlation where function and the sampling frequency, so it can be treated as a known quantity after the estimation of the correlation function. Particularly, for PSK signals with code-rate of , can be computed directly as follows:

If the ambiguity occurs also for the incident signals holistically, then the following equality can be obtained by dividing (24) into the higher and lower spectral parts,

(33)

(29) Combining (28) and (29) yields

(30) As the array geometry is designed to avoid ambiguity within holds for all the frequency band , . Thus, one can easily demonstrate that (30) is never satisfied following the derivation of Theorem 1. Therefore, any potential ambiguity within the lower band is got rid of due to the nonambiguity guarantee within the higher band. V. COMPUTATION OF THE PERTURBATION VARIANCES The fitting error thresholds in (16) and (20) depend on the perturbation levels of the measurement vectors, so we present a method for computing their variances in this section to facilitate the implementation of W-CMSR. The perturbation variance expectations (or upper bound) given in (15) and (19) can be rewritten as

(31) and

In addition, the expression of in (2) indicates that , thus can be taken as an estimate , denote it by . Furthermore, the noise variof ance estimate is also required for the computation of (31) and (32). In the narrowband setting, the average of the smaller eigenvalues of the array output covariance matrix corresponding to the noise subspace provides a good estimate for noise variance. But for wideband signals, the signal energy disperses to more subspaces than the signal number, so we eliminate all those larger eigenvalues (generally more than true signal number) and take the average of the smallest ones as the noise variance estimate, i.e., (34) where is the largest eigen-value of the array output cois the signal number estimate derived variance matrix, and from general narrowband source enumeration methods, such as MDL [22]. and into (31) and (32) yields the perSubstituting , turbation variance expectation (or upper bound) of the measurement vectors in GLA and ULA. After that, adequate weighting factors of and should be chosen to obtain the fitting error thresholds in (16) and (20). Finally, (11) and (18) can be solved for wideband DOA estimation. VI. SIMULATION RESULTS In this section, we carry out simulations to analyze the performance of W-CMSR in wideband DOA estimation, and also compare it with that of previous methods, including CSSM, Capon, L1-SVD,2 and JLZA-DOA.3 Two arrays with the same aperture of six half-wavelengths are used in the simulations, one is the seven-element ULA and the other is the four-element MRA, whose elements are located at 0, 1, 4, and 6 spacing 2The 3The

source code of L1-SVD is provided by Malioutov [7]. source code of JLZA-DOA is provided by Hyder [8].

4262

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

Fig. 1. Spectra on 7-ULA interspaced by half-wavelength with respect to the highest signal frequency. (a) Capon; (b) L1-SVD; (c) JLZA-DOA; (d) CSSM, W-CMSR and spectrally averaged (a), (b), and (c).

units [14]. If not specially stated, the half-wavelength is set with respect to the highest signal frequency to avoid ambiguity in all the methods. For notation simplification, we denote those two arrays by 4-MRA and 7-ULA, respectively. When 4-MRA is used, no repeated off-diagonal elements exist in the perturbation-free covariance matrix, so we choose the method of W-CMSR(GLA) for DOA estimation. But in 7-ULA, both methods of W-CMSR(GLA) and W-CMSR(ULA) are usable. Among the methods considered, W-CMSR, Capon and JLZA-DOA do not require the a priori information of source number, but in CSSM and L1-SVD the source number is assumed known beforehand. The optimization tool of SeDuMi [23] is introduced to solve the L1-norm-based optimization problems of (11), (18) and also the one in L1-SVD. In all of spatial scope is sampled with an the methods, the interval of 1 . The snapshots are divided into sections of 64 snapshots each in CSSM, Capon, L1-SVD, and JLZA-DOA for temporal-spectral transformation and spectral snapshot collection. and When establishing the overcomplete dictionaries of for W-CMSR, (9) is used for correlation computation when the incident signals are PSK ones with known code-rate, while (10) is used otherwise. And when computing the correlations with (10), we first estimate the signal power spectrum from the array output with the function periodogram in Matlab, and then use IFFT to compute the correlation at certain timedelays. In the simulations, we only use the array output at the to save computation load. first sensor to estimate

In Part A of this section, we first give an intuitionistic spectrum comparison between W-CMSR and previous methods. Then we relax the inter-spacing of the array from half-wavelength with respect to the highest signal frequency to that with respect to the lowest signal frequency, so as to demonstrate the half-wavelength relaxation property of W-CMSR, and also provide guidelines for array design in subsequent simulations. The performances of those methods are then compared in at different the statistical sense in Part B. We set and values for W-CMSR(GLA) and W-CMSR(ULA) in Part B, and carry out simulations for two signals at various angle separations to test the impacts of array geometry, fitting error threshold and prior spectrum information on the performance of W-CMSR. Empirical values of and can be derived from this group of simulations. After that, W-CMSR is compared with other methods under different inter-spacings and various snapshots and SNR, and comparison is also made with respect to computational efficiency. Finally, the ability of W-CMSR in separating multiple signals is demonstrated in Part C. A. Spectra and Half-Wavelength Relaxation Suppose that two 0-dB BPSK signals with central frequency of 70 MHz and bandwidth of 40% impinge onto a 7-ULA from the directions of 30 and 40 , 64 snapshots are collected. The methods of CSSM, Capon, L1-SVD, JLZA-DOA, and are used to estimate the W-CMSR(ULA) with directions of those signals, the signal spectrum is assumed unknown for W-CMSR. First, fix the inter-spacing of the array

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

4263

Fig. 2. Spectra on 7-ULA inter-spaced by half-wavelength with respect to the lowest signal frequency. (a) Capon; (b) L1-SVD; (c) JLZA-DOA; (d) CSSM, W-CMSR and spectrally averaged (a), (b), and (c).

at half-wavelength with respect to the highest signal frequency, the spectra of the five methods in a randomly chosen trial are given in Fig. 1. As Capon, L1-SVD and JLZA-DOA decompose the incident signals into narrowband components first, and obtain spatial spectra within each frequency bin, we retain the whole spectral-spatial spectra of those three methods in Fig. 1(a)–(c) to give a more compositive view of the possible shortcomings of these methods. Those 3-D spectra indicate that the peak locations and amplitudes within different frequency bins are generally quite different, thus no clear-cut guideline is available for deriving integrative wideband DOA estimates from them. We then average the 3-D spectra of Capon, L1-SVD and JLZA-DOA in the spectral domain and reprint the corresponding 2-D spectra in Fig. 1(d) for a more convenient comparison with CSSM and W-CMSR. The results in Fig. 1 indicate that CSSM, Capon, L1-SVD and JLZA-DOA all fail to separate the two signals, and an amount of pseudo-peaks exist in the spectrum of L1-SVD. On the contrary, W-CMSR gains two significant peaks at 31 and 37 , with a maximal direction estimation bias of 3 . Then, we enlarge the inter-spacing of the array to half-wavelength with respect to the lowest signal frequency, which is one and a half times of that in the above simulation, and fix the other settings unchanged. The spectra of the five methods in a randomly chosen trial are given in Fig. 2. Among those subfigures, (a)-(c) illustrate the 3-D spectra of Capon, L1-SVD, and JLZA-DOA, respectively, and they are averaged spectrally

and reprinted in Fig. 2(d) to make a more direct comparison with CSSM and W-CMSR. As the conventional half-wavelength constraint is not satisfied, significant ambiguity occurs in the spectra of Capon and L1-SVD, together with an amount of pseudo-peaks in the spectrum of L1-SVD. Ambiguity also emerges near -90 in CSSM, as is shown in Fig. 2(d). In addition, all of those three methods fail to separate the two incident signals. Due to the synthesized exploitation of the signal energy within the band, no ambiguity occurs in JLZA-DOA or W-CMSR. However, JLZA-DOA still fails to separate the two signals, and only W-CMSR succeeds again. The two spectrum peaks of W-CMSR are located at 28 and 39 , respectively, with a maximal direction estimation bias of 2 , which is smaller than the 3 bias in the previous simulation. B. Statistical Performance Comparison This part compares the performance of W-CMSR and several existing DOA estimation methods, including CSSM, Capon, L1-SVD and JLZA-DOA, in the statistical sense for signals of various separation angles, snapshot numbers and signal-to-noise ratios (SNRs). The computational efficiencies of those methods at different snapshot numbers and signal bandwidths are also compared. In the first group of simulations, we assume that two 0-dB BPSK signals with bandwidth of 20% impinge onto the 4-MRA and 7-ULA, 256 snapshots are collected, and the direction of the first signal is fixed at 10 , while that of the second signal

4264

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

Fig. 3. Separation probabilities of W-CMSR, together with CSSM, Capon, L1-SVD, and JLZA-DOA, for two BPSK signals at various angular separations with DOA1 = 10 , SNR = 0 dB, = 256 and 20% bandwidth. (a) W-CMSR(GLA) on 4-MRA, code-rate unknown; (b) W-CMSR(GLA) on 4-MRA, code-rate known; (c) W-CMSR(GLA) on 7-ULA, code-rate unknown; (d) W-CMSR(GLA) on 7-ULA, code-rate known; (e) W-CMSR(ULA) on 7-ULA, code-rate unknown; and (f) W-CMSR(ULA) on 7-ULA, code-rate known.

N

varies from 17 to 30 . Different implementations of W-CMSR with various weighting factors are carried out, together with CSSM, Capon, L1-SVD, and JLZA-DOA. For each angle separation of the two signals, the separation probabilities of all the methods are derived from 1000 trials. Successful separation of CSSM and W-CMSR is defined if the DOA estimation biases of both signals do not exceed 3 , and the amplitude of the highest pseudo-peak does not exceed half of that of the lowest signal peak. But for Capon, L1-SVD, and JLZA-DOA, as

they obtain spectra in discrete frequency bins, and the results in Fig. 1(a)–(c), and Fig. 2(a)–(c) show that the locations and amplitudes of the spectrum peaks within different frequency bins are generally quite dissimilar, thus we carry out the (less than 3 bias and lower than -3 dB pseudo-peaks) criterion within each bin for separation evaluation since there exists no clear-cut guidelines for integrating those narrowband spectra. If the two incident signals are separated in more than half of the frequency bins, we take it for successful separation.

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

4265

Fig. 4. Separation probabilities of CSSM, Capon, L1-SVD, JLZA-DOA, and W-CMSR for two BPSK signals with 20% bandwidth on 7-ULA, signal directions , SNR varies from -5 dB to 10 dB. are fixed at 10 and 20 . (a) SNR = 0 dB, snapshot number varies from 64 to 1024; (b)

N = 256

Firstly, W-CMSR(GLA) is carried out on 4-MRA with weighting factors set at 0.5, 1, 1.5, 2, and 2.5. The over-complete dictionary is established according to (10) without any prior spectrum information or to (9) with known code-rate. The separation probabilities of W-CMSR at different angle separations are shown in Fig. 3(a) for unknown code-rate and Fig. 3(b) for known code-rate, respectively. Then keep the other settings unchanged and use a 7-ULA for DOA estimation, both implementations of W-CMSR(GLA) and W-CMSR(ULA) are carried out. The weighting factors of W-CMSR(GLA) are chosen the same as those in 4-MRA, the corresponding separation probabilities are given in Fig. 3(c) and (d) when the signal code-rate is assumed unknown and known, respectively. The weighting factors of W-CMSR(ULA) are set at 0.3, 0.5, 0.7, 0.9, and 1.1, and its separation probabilities based on (10) and (9) are plotted in Fig. 3(e) and (f), respectively. The dashed curves in Fig. 3(a)–(f) corresponds to the separation probabilities of CSSM, Capon, L1-SVD, and JLZA-DOA. It can be concluded from Fig. 3 that when the fitting error threshold of W-CMSR is set within a certain scope, enhanced superresolution can be obtained from lower thresholds, but false detection of different degrees emerges for even well-separated signals when the threshold is set too low. If higher threshold is chosen, superresolution is weakened, but it is more potential to obtain 100% separation probability for well-separated signals. When the code-rate of the incident BPSK signals is not known beforehand, well-tuned W-CMSR(GLA) surpasses the other methods on 4-MRA, and such superiority is much more significant for W-CMSR(ULA) on 7-ULA. However, W-CMSR(GLA) does not perform well in 7-ULA in most of the given settings. A heuristic interpretation for such performance deterioration from W-CMSR(ULA) to W-CMSR(GLA) on ULA is that the more repeated elements in makes the distribution of flat, and it becomes difficult to give attention to both better superresolution and weaker pseudo-peaks simultaneously. Therefore, we recommend combining the repeated elements in if there are any. The comparison between the three pairs of counterparts when the code-rate is assumed known and unknown indicates that, the exploitation of the prior information of signal spectrum helps to largely improve the performance of W-CMSR. With the signal code-rate assumed known, W-CMSR(GLA)

and W-CMSR(ULA) surpasses the other four methods largely on 4-MRA and 7-ULA, respectively, and the 7-ULA-based W-CMSR(GLA) also gains some superiority over the other methods if the threshold is appropriately set. Among the other four methods, L1-SVD performs best in superresolution, but as there are too many pseudo-peaks in its spectrum, it can not separate well-separated signals with 100% probability. Coherent subspace-based CSSM exceeds incoherent Capon in most of the settings, while JLZA-DOA does not perform well all along. We then fix the directions of the two BPSK signals with 20% bandwidth at 10 and 20 , and test the adaptability of W-CMSR to snapshot number and SNR. The 7-ULA is used for DOA estimation, and W-CMSR is implemented via W-CMSR(ULA) with . No prior spectrum information is introduced for W-CMSR unless otherwise stated. First fix the SNR of both signals at 0 dB and vary the snapshot number from 64 to 1024, CSSM, Capon, L1-SVD, JLZA-DOA, and W-CMSR are carried out for comparison, and their probabilities of separating those two signals are shown in Fig. 4(a). According to the simulation results in Part A of this section, only JLZA-DOA and W-CMSR among the five methods survive from potential spatial ambiguity when the inter-spacing is relaxed to half-wavelength with respect to the lowest signal frequency, so we perform them in such settings as well to show the impact of array enlargement on them. The variations of JLZA-DOA and W-CMSR on enlarged array are labeled as JLZA-DOA(AE) and W-CMSR(AE) (AE: Array Enlarged). The impact of prior spectrum information on W-CMSR is also tested, and the variation method is labeled as W-CMSR(SK) (SK: Spectrum Known). Then we fix the snapshot number at 256 and vary the SNR of both signals from -5 dB to 10 dB, the probabilities of those five methods and the variations of JLZA-DOA and W-CMSR in separating the two signals are shown in Fig. 4(b). Those two groups of simulation results indicate that, W-CMSR is much more adaptable to fewer snapshots and lower SNR than the other four methods, and enlarged array aperture and prior spectrum information help to further improve its performance. Among the other four methods, L1-SVD performs best and CSSM exceeds Capon and JLZA-DOA in most of the settings, JLZA-DOA and its variation on enlarged array both fail to perform well even when the snapshots are sufficient or the SNR is high. In-depth analysis of the JLZA-DOA spectrum shows that,

4266

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

TABLE I COMPUTATION TIME COMPARISON

Fig. 5. Separation of six simultaneous wideband signals with W-CMSR. (a) 4-MRA and (b) 7-ULA.

but linearly with signal bandwidth, it is computationally most demanding when the snapshot number is 64 or 256. The computational efficiency of W-CMSR is not affected by the signal bandwidth, and deteriorate slightly when the snapshot number is too large due to the periodogram and IFFT processes during the establishment of the over-complete dictionary. W-CMSR is faster than L1-SVD in all the settings considered, and is also faster than JLZA-DOA when the bandwidth is large or the snapshots are sufficient. In addition, the signal spectrum is assumed unknown beforehand for W-CMSR in this simulation, and we compute the over-complete dictionary online in each trial, this process greatly aggravates the computation load. We have concluded from previous simulations in this part that the establishment of the dictionary takes more than 80% of the computation time in W-CMSR. Therefore, if prior spectrum information is available, W-CMSR will outperform L1-SVD and JLZA-DOA significantly in computational efficiency. Furthermore, it is more important to note that W-CMSR performs much better than the other four methods in DOA estimation. C. Separation of Multiple Signals

the true spectrum peaks often fork at the top, which may puzzle the user in judging whether there is one signal or two adjacent signals (see Fig. 5 for details). This phenomenon has blocked JLZA-DOA from gaining satisfying DOA estimation performance, and enlarged array aperture does not help to eliminate, or even weaken, the influence of the forked peaks. Finally, we compare the average computation time of those five methods for various snapshot numbers and signal bandwidths. Two 0 dB signals with directions of 10 and 20 are assumed to impinge simultaneously, and the incident signals are formed as mixtures of equal-power sinusoids within the signal band, instead of BPSK ones in the other simulations, to make the bandwidth more conveniently adjustable. The signal bandwidth is set at 20%, 40%, and 80%, and the snapshot number is set at 64, 256, and 1024. For each bandwidth-snapshot number pair, 1000 trials are carried out to obtain the average computation time of those five methods, and the results are given in Table I. The results show that CSSM and Capon are the two fastest methods. Among the three sparsity-inducing methods, JLZA-DOA is faster at narrower bandwidth and fewer snapshots, but its computation time increases linearly with bandwidth and snapshot number. When the snapshot number is 1024, JLZA-DOA needs up to 18.56 s to perform direction estimation for two signals with 80% bandwidth once! The computation time of L1-SVD increases slightly with snapshot number,

In this part, we carry out simulations to demonstrate the ability of W-CMSR in separating multiple wideband signals. As in the above simulations, 4-MRA and 7-ULA are used. According to Theorem 1, both arrays can separate up to six simultaneous signals. So we suppose that six 5 dB BPSK signals with 20% bandwidth impinge onto the arrays from -55 , -25 , 3 , 10 , 35 , and 60 , respectively, with a minimum angle separation of 7 , 512 snapshots are collected. The code-rate of the incident signals is assumed unknown and the over-complete dictionary of W-CMSR is established according to (10). In 4-MRA, the method of W-CMSR(GLA) is used with , and the spectra of Capon, L1-SVD, JLZA-DOA and W-CMSR are shown in Fig. 5(a), in which the spectra of Capon, L1-SVD and JLZA-DOA have been averaged spectrally as that in Figs. 1(d) and 2(d). As the incident signal number is larger than the sensor number, CSSM is excluded since there is no noise subspace any more. In 7-ULA, the method of W-CMSR(ULA) is used with , and the (spectrally averaged if necessary) spectra of those four methods, together with CSSM, are shown in Fig. 5(b). The simulation results indicate that W-CMSR can successfully separate those six simultaneous signals on both 4-MRA and 7-ULA, while the other four methods all fail, especially in separating the 3 and 10 signals. It can also be seen that the JLZA-DOA spectrum peaks fork at the top for even the far-apart signals. This phenomenon illustrates why JLZA-DOA

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

did not perform well in most of the settings in Part B. Further simulations show that seven signals can never be separated by any of those methods on both arrays. This conclusion well supports Theorem 1 in Section IV, and it is much inspiring that W-CMSR succeeds to separate more signals than sensors with well-designed arrays, such as MRA. As far as we know, none of the previous wideband DOA estimation methods has achieved such watershed, no matter subspace-based ones [1]–[5] or sparsity-based ones [6]–[8].

4267

(35) Denote

, then

VII. CONCLUSION A wideband DOA estimation method, named W-CMSR, is proposed in this paper. This method does not need the a priori information of the incident signal number, and does not introduce any spectral decomposition or focusing procedures, so it is more adaptable to noncooperative applications and never run into the trouble of multi-subband fusion to get final DOA estimates. Simulation results show that W-CMSR surpasses CSSM, Capon, L1-SVD, and JLZA-DOA in DOA estimation if the fitting error threshold is set appropriately. Enlarged array aperture and prior signal spectrum information further improves the performance of W-CMSR. Theoretical analysis and simulation results discovered that W-CMSR has an intrinsic property to make good use of the array geometry, such as MRA, to separate even more simultaneous signals than sensors. It is also demonstrated that the half-wavelength constraint for nonambiguous DOA estimation is relaxed from the highest to the lowest signal frequency in W-CMSR. APPENDIX A PERTURBATION VARIANCE OF THE MEASUREMENT VECTOR IN GENERAL ARRAYS element of the covariance matrix estimated from The snapshots is

(36) Define , then the expectation of the perturbation variance of this element is (37) Combining the zero-mean and Gaussian distribution of the incident signals and additive noise, it can be concluded that

(38) And the first nonzero entity of in (37) can be computed as (39) at the bottom of the page. As , we can make the following approximation: (40)

(39)

4268

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

where depends on the signal correlation function and the sampling frequency. Thus, (39) is simplified to

In addition,

(41)

(46)

(48)

(49)

(50)

(51)

LIU et al.: DoA ESTIMATION OF WIDEBAND SIGNALS VIA COVARIANCE MATRIX SPARSE REPRESENTATION

4269

APPENDIX B PERTURBATION VARIANCE OF THE MEASUREMENT VECTOR IN ULA

(42)

yields Substituting (36) into See (46) at the bottom of the previous page. , then it can be concluded Define from the properties of zero-mean and Gaussian distribution of the incident signals and additive noise that (47) And the other four nonzero elements can be computed as (48)–(51) at the bottom of the previous page. Thus

Similar derivations can be carried out to conclude that

(43) (44) Synthesizing (37)–(44) yields

Equation (45) indicates that element indexes (15).

and

(45)

(52)

is independent of the

So the perturbation variance of the measurement vector satisfies (53) at the bottom of the page. The approximation in the fourth row of (53) is introduced because for

, so substituting (45) into (14) yields

(53)

4270

elements in the set equal , while the others distribute symmetrically on both , and the correlation function of most practical sigsides of . The approximation nals decreases on average besides in the fifth row bases on (40). ACKNOWLEDGMENT The authors would like to thank the anonymous reviewers for valuable comments and suggestions that have helped to vastly improve the contents and presentations of this paper. The authors are grateful to D. Malioutov and M. M. Hyder for providing the source codes of L1-SVD and JLZA-DOA, and appreciate D. Malioutov also for the fruitful discussions on L1-normbased methods. REFERENCES [1] M. Wax and T. Kailath, “Spatio-temporal spectral analysis by eigenstructure methods,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-32, no. 4, pp. 817–827, Aug. 1984. [2] H. Wang and M. Kaveh, “Coherent signal-subspace processing for the detection and estimation of angles of arrival of multiple wideband sources,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-33, no. 4, pp. 823–831, Aug. 1985. [3] M. A. Doron and A. J. Weiss, “On focusing matrices for wideband array processing,” IEEE Trans. Signal Process., vol. 40, no. 6, pp. 1295–1302, Jun. 1992. [4] Y. S. Yoon, L. M. Kaplan, and J. H. McClellan, “TOPS: New DOA estimator for wideband signals,” IEEE Trans. Signal Process., vol. 54, no. 6, pp. 1977–1988, Jun. 2006. [5] E. D. D. Claudio and R. Parisi, “WAVES: Weighted average of signal subspaces for robust wideband direction finding,” IEEE Trans. Signal Process., vol. 49, no. 10, pp. 2179–2190, Oct. 2001. [6] J. J. Fuchs, “On the application of the global matched filter to DOA estimation with uniform circular arrays,” IEEE Trans. Signal Process., vol. 49, no. 4, pp. 702–709, Apr. 2001. [7] D. Malioutov, M. Cetin, and A. S. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 3010–3022, Aug. 2005. [8] M. M. Hyder and K. Mahata, “Direction-of-arrival estimation using a mixed L2,0 norm approximation,” IEEE Trans. Signal Process., vol. 58, no. 9, pp. 4646–4655, Sep. 2010. [9] B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM J. Comput., vol. 24, pp. 227–234, 1995. [10] J. A. Tropp, “Just relax: Convex programming methods for identifying sparse signals in noise,” IEEE Trans. Inf. Theory, vol. 52, no. 3, Mar. 2006. [11] J. A. Tropp and S. J. Wright, “Computational methods for sparse solution of linear inverse problems,” Proc. IEEE, vol. 98, no. 6, pp. 948–958, Jun. 2010. [12] M. M. Hyder and K. Mahata, “A robust algorithm for joint-sparse recovery,” IEEE Signal Process. Lett., vol. 16, no. 12, pp. 1091–1094, Dec. 2009. [13] Z. M. Liu, Z. T. Huang, and Y. Y. Zhou, “Source number detection and direction estimation via sparsity-inducing representation of the array covariance matrix,” IEEE Trans. Aerosp. Electron. Syst., to be published. [14] A. T. Moffet, “Minimum-redundancy linear arrays,” IEEE Trans. Antennas Propag., vol. AP-16, no. 2, pp. 172–175, Mar. 1968. [15] J. G. Proakis, Digital Communications, 4th ed. New York: McGrawHill, 2001.

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 59, NO. 9, SEPTEMBER 2011

[16] H. Krim and M. Viberg, “Two decades of array signal processing research: The parametric approach,” IEEE Signal Process. Mag., vol. 13, no. 4, pp. 67–94, Jul. 1996. [17] B. P. Lathi, Signal Processing and Linear Systems. London, U.K.: Oxford Univ. Press, 1998. [18] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Rev., vol. 43, no. 1, pp. 129–159, 2001. [19] M. Zibulevsky and M. Elad, “L1–L2 optimization in signal and image processing,” IEEE Signal Process. Mag., pp. 76–88, May 2010. [20] D. L. Donoho, S. Mallat, R. Sachs, and Y. Samuelides, “Locally stationary covariance and signal estimation with macrotiles,” IEEE Trans. Signal Process., vol. 51, no. 3, pp. 614–627, Mar. 2003. [21] D. L. Donoho and J. Tanner, “Sparse nonnegative solution of underdetermined linear equations by linear programming,” Proc. Nat. Acad. Sci., vol. 102, no. 27, pp. 9446–9451, Jul. 2005. [22] P. Stoica and Y. Selen, “Model-order selection: A review of information criterion rules,” IEEE Signal Process. Mag., vol. 21, no. 4, pp. 36–47, Jul. 2004. [23] J. S. Sturm, Using SeDuMi 1.02, A Matlab Toolbox for Optimization Over Symmetric Cones. Tilburg, The Netherlands, Dept. Econometrics, Tiburg Univ., 2010 [Online]. Available: http://fewcal.kub.nl/ ~strum

Zhang-Meng Liu (S’09) received the Bachelor’s degree in electronic science and engineering from the National University of Defense Technology, Changsha, Hunan, China, in 2006. Currently, he is working towards the Ph.D. degree in the National University of Defense Technology, Changsha, Hunan, China. His research interests include array signal processing in radar and communication applications, sparse reconstruction.

Zhi-Tao Huang received the B.S. and Ph.D. degrees in information and communication engineering from the National University of Defense Technology, Changsha, Hunan, China, in 1998 and 2003, respectively. He is currently a Professor with the School of Electronic Science and Engineering, National University of Defense Technology, Changsha, Hunan, China. His research interests include radar and communication signal processing, array signal processing.

Yi-Yu Zhou received the B.S., M.S., and Ph.D. degrees in information and communication engineering from the National University of Defense Technology, Changsha, Hunan, China, in 1982, 1985, and 1992, respectively. He is currently a Professor with the School of Electronic Science and Engineering, National University of Defense Technology, Changsha, Hunan, China. His research interests include radar signal processing and passive location.