DOA Estimation Of Wideband Sources Using A ... - Semantic Scholar
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
619
DOA Estimation of Wideband Sources Using a Harmonic Source Model and Uniform Linear Array Monika Agrawal and Surendra Prasad, Senior Member, IEEE Abstract— We consider here the problem of estimation of DOA’s of multiple wideband sources incident on a uniform linear array (ULA) in the presence of spatially and temporally white Gaussian noise (WGN). The approach presented here builds up on the IQML algorithm suggested by Bresler and Macovski for the case of narrowband DOA estimation. It is shown that the concept of ARMA model for the observed data vector for the narrowband case can be generalized to model an appropriately stacked, space-time data vector obtained by combining the spacetime samples. The coefficients of the corresponding 2-D predictor polynomial can be used to represent the null subspace of the wideband array steering matrix, and rooting of the polynomial at each frequency, separately, gives the DOA estimates. These separate estimates at multiple frequencies are combined into a single DOA estimate in a least squares sense. This leads to the formulation of an IQML-like procedure for the spatial parameter estimation of wideband sources. Like its narrowband counterpart, the proposed approach is applicable to both noncoherent and coherent sources. The performance of the proposed method is studied via extensive computer simulations and by comparison with the CR bounds. Index Terms— Array signal processing, broadband sources, direction-of-arrival estimation, maximum likelihood techniques, special ARMA model, uniform linear array.
I. INTRODUCTION
T
HE PROBLEM of estimating the source parameters by processing data from a sensor array has received much attention in recent years. In so far as the estimation of the parameters of the narrowband sources is concerned, the theory is well established, and a large body of literature exists. The subspace-based methods have evolved into a powerful class of methods for direction-of-arrival (DOA) estimation of narrowband sources. These methods have also helped the formulation and solution of the maximum likelihood (ML) approach to the problem. For wideband sources, however, the literature is somewhat scanty, although a number of researchers have focussed their attention to this problem in recent years. Wax et al. [1], Wang and Kaveh [2], Buckley and Griffiths [3], and Morf and Su [4] are among the earlier researchers in this field who suggested some useful and interesting solutions for the problem of DOA estimation of wideband sources. Both incoherent aggregation suggested in [1] and coherent signal subspace (CSS) technique proposed in [2] represent in Manuscript received November 6, 1997; revised August 21, 1998. The associate editor coordinating the review of this paper and approving it for publication was Dr. Eric Moulines. The authors are with the Department of Electrical Engineering, Indian Institute of Technology, Delhi, New Delhi, India. Publisher Item Identifier S 1053-587X(99)01328-8.
some way adaption of narrowband processing techniques to the wideband case. Bienvenu [5] and Buckley and Griffiths [3] have shown how to divide the spatio-temporal observation space into signal and noise subspaces and have suggested search functions to estimate the DOA’s. Grenier [6] has generalized this notion to a so-called filter-spatial domain. Morf and Su [4] have suggested a modal decomposition of the signal subspace for estimating the DOA’s. In this paper, we are concerned with the solution of the wideband DOA estimation problem based on maximum likelihood principles. The problem of maximum likelihood estimation has, of course, been addressed earlier by Doron et al. [6], who have proposed several search functions to obtain the ML estimate of DOA for wideband sources. However, such an exhaustive search is computationally expensive and impractical for real-time applications. In this respect, the approach proposed by Bresler and Macovski [8] for the narrowband case, and a uniform linear array (ULA), is of particular interest, as it formulates an iterative method called the iterative quadratic maximum likelihood (IQML) to obtain the ML solution. The IQML method involves the creation of a Toeplitz matrix that spans the subspace orthogonal to the signal subspace and can be parameterized by the coefficients of a prediction polynomial, whose roots yield the DOA estimates. In effect, therefore, use of this parameterized model converts the problem of estimating the DOA’s to that of estimating the coefficients of this model from which DOA’s can be easily calculated. Clark and Scharf [9] have extended this work to twodimensional (2-D) damped harmonic signals in additive Gaussian noise which, in principle, can take care of the wideband case. However, as pointed out by them, it is not possible to directly generalize the one-dimensional (1-D) results to the 2-D case due to the lack of a fundamental theorem of algebra for 2-D polynomials. They suggested a solution to the problem by using a prediction polynomial in one of the two dimensions and a polynomial that interpolates between the two. This paper builds on the ideas available for the parameter estimation of multiple narrowband sources in WGN to develop pseudo-maximum likelihood parameter estimators for the wideband case by adapting the harmonic model of Kay et al. [10] for wideband sources. It is shown that the concept of an ARMA model for the observed data vector used in the narrowband case [8] can be generalized to model an appropriately stacked, space-time data vector obtained by combining the spatio-temporal samples. As in the narrowband case, it is shown that the coefficients of the resulting ARMA model can be used to represent the null subspace of the wideband
1053–587X/99$10.00 1999 IEEE
620
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
array steering matrix associated with the model of the stacked observations for the harmonic signal model considered here. The resulting 2-D predictor polynomial is rooted for each frequency separately to give the DOA estimates. This leads us to the formulation of an IQML-like procedure for the parameter estimation of wideband sources that is similar to that for the narrowband case suggested in [8]. Thus, the difficulties associated with the 2-D model analysis due to lack of a theorem of algebra for 2-D polynomials or with pairing the roots when the 2-D problem is addressed by breaking it into two 1-D problems (as in [11]), are avoided in this approach.1 The solution, however, falls short of providing a single ML estimate for each source. This is because the IQML search is conducted over an expanded signal space, which in turn leads to generation of separate DOA estimates for different frequencies. For applications where it would be more desirable to obtain a single DOA estimate for each wideband source, we suggest a modification of the IQML search whereby the separate estimates at multiple frequencies can be combined into a single estimate of the DOA in an appropriately defined least squares sense. It is seen that the procedure leads to very good estimates, as evidenced by comparison with the Cram´er–Rao bound. This paper has been organized as follows. In Section II, we formulate the problem and develop its mathematical model. In Section III, we develop the ML criterion in a parametric framework. Sections IV and V constitute the main body of our work, where we unfold our solution. Results and simulation studies are presented in Section VI. II. MODELING
AND
Let be the total number of wideband sources present. In this paper, we make the usual assumption of to be known. may be estimated using the methods given Alternatively, in [2]. Thus, at the th sensor, the th time sample of the output signal corresponding to th snapshot is given by (2) is spatially and temporally white Gaussian where such that noise of variance (3) Define
(4) where
is the number of sensors present in the ULA. Let (5)
be the temporal and spatial frequency variables, respectively. Using this notation and combining (1)–(3), we can write
(6) The above equation can be rewritten in a vector form as
PROBLEM FORMULATION
As mentioned, we model a broadband signal here as the sum of a finite number of monochromatic propagating plane waves. Each source is assumed to have a spectral support of Hence, the th time sample of the signal from the th (relative to the broad side source impinging from direction of the array) at the th sensor of a ULA corresponding to the th snapshot can be expressed as [10]
(7) where
(8) and
(1) where ’s are harmonically related frequencies representing is the the spectral components of the sources. In addition, Fourier coefficient of th signal corresponding to , the th , where frequency, array element spacing; velocity of propagation of plane waves; number of monochromatic waves used to model the broadband source. 1 It has been brought to our attention by one of the reviewers that there is some similarity of our approach to that of N. Chotikakamthorn, as reported in his doctoral work at the Department of Electrical Engineering, Imperial College of Science, Technology and Medicine, University of London [22]. However, at the time of submitting this paper, we did not have access to this dissertation; we were, therefore, unable to comment on the precise relationship of our work with that of [22].
(9) to be a block It is convenient to regard the vector may Vandermonde vector. A block Vandermonde matrix now be defined as (9a) , and are where the temporal and spatial parameters of the signal. Defining , the system equation (7) can be rewritten in a concise form as (10) ’s and In this formulation, the unknown parameters are ’s. Hence, our estimation problem should be formulated as follows.
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
Given snapshots of the data vector , estimate ’s and ’s, where ’s are embedded in the steering matrix , and ’s are the components of the vector However, this estimation would require the solution of a highly nonlinear parameters optimization problem in terms of the Alternatively, we can consider the problem of esticontains the mating the spatial vector , whose component required DOA information. Even here, ideally, the estimation problem should be posed in such a way that the components of the estimated are related to each other via (5), i.e., In the developments that follow in Sections III and IV, it is seen that there is no obvious algebraic procedure for imbedding such a constraint on the parameter vector It is convenient, therefore, to split this problem into two steps. In the first step, we consider the modified problem of independently at all relevant frequencies, i.e., estimating Using (5), this step would lead to for separate estimates of : one corresponding to each frequency in the harmonic model. In the second step, we exploit the implied structural relationship between by (5) to obtain a single estimate for DOA of each wideband source in a least squares sense. These two steps are considered in more detail in the sequel. III. MAXIMUM LIKELIHOOD CRITERION As discussed in the proceeding paragraph, the first step of the problem is concerned with estimation of the parameter vector More precisely, given snapshots of the data vector , estimate ’s and ’s, where ’s are embedded in the ’s are the components of the vector steering matrix , and The solution of this problem is taken up in this section. It has been shown that under WGN, ML estimators and least square estimators are equivalent [12]. Hence, ML estimate of the signal parameters can be obtained by solving the nonlinear least squares problem (11) is the Euclidean norm. This expression may be where further simplified by first substituting the least square estimate in terms of as of (12) The least square estimate of the can then be shown to be given by [13]
where vector
arg
arg
tr
(13)
arg
tr
(14)
or, equivalently arg
and are the projection matrices for projection Here, onto the column space of and onto its orthogonal complement, respectively, and are given by
(15)
where tr is the trace operator, and given by covariance matrix of
621
is the estimated
(16) being the number of snapshots of the observation vector Having obtained the extended spatial parameter vector (of ) by solving (13) or (14), our next step is to look dimension for a -dimensional vector of unknown directions that best fits the estimated vector in the least squares sense. In order to do this, we recall that for the th source, the components of the vector are given by Assuming that the ML estimation of would bring about a decoupling of the estimates of its subvectors from each other, a reasonable (although, in principle, not exact) formulation of the least squares problem to estimate can be stated as follows. For the estimated such that vector , find an angle (17) is minimum. Such a least squares fit helps to combine the independently obtained DOA estimates at the various frequencies into a single effective estimate for each source. In the next section, we consider first the issue of estimating the spatial vector via a characteristic polynomial represenThe tation of the broadband, spatio-temporal observations detailed algorithm that combines the two steps of estimating and then in an iterative search framework is taken up in Section V. IV. ML CRITERION AND THE MODEL PARAMETERS It has been shown in the literature for the narrowband case [8] that the output vector observed on a ULA in the presence of WGN obeys a special ARMA model and that the ML estimates of its parameters are directly related to those of its coefficients. In effect, use of this model converts the problem of estimating the DOA’s to that of estimating the special ARMA parameters from which the DOA’s can be easily derived. Taking a clue from this equivalence, we consider the possibility of formulating a similar, equivalent for the ARMA model for the spatio-temporal samples broadband case, as given by (6). In view of the spatio-temporal , it is proposed that we use an nature of the samples expanded form of the special ARMA model for these samples. More specifically, we establish the following results for the space-time vector Theorem 1: The space-time vector admits the ARMA representation given by (18) where the set of coefficients ’s are assumed to form a 2-D linear predictor polynomial given by (18a)
622
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
with and denoting the temporal and spatial complex frequency variables, respectively. Proof: Appendix A outlines a proof of this equivalent by adopting the methodology of Ulrych and description of Clayton [14]. The polynomial may be regarded as the linear predictor polynomial for the noiseless signal component of (2). This model can be considered to be a generalization of the Ulrych and Clayton model [14] for the case of 2-D harmonic of series. For convenience, we define the vectors and , respectively, as lengths
where each block is a circulant matrix of size
given as
(20a)
Then,
can be written as (21) be a column of and
Proof: Let to see that for
Then, it is easy
(19a) and
(22)
(19b) Here, elements of the vector polynomial for each as
form a 1-D linear predictor
(19c) is the th element of the vector It follows where ’s of this polynomial are known, that if the coefficients ’s can be estimated by rooting the spatial frequencies This result forms the basis of our proposed approach for the parameter estimation problem for wideband (more precisely multichromatic) plane waves in spatially and temporally white Gaussian noise. The method is similar to that for the narrowband case and consists of following steps. A) Establish a relation to represent the null space of the array steering matrix associated with the model of the space-time stacked vector , as defined by (10), in terms of the coefficients of the special ARMA model. B) Formulate the minimization problem in terms of a constrained quadratic minimization procedure for obtaining the ARMA coefficients Like the narrowband case, once the coefficients of the are estimated, the unknown DOA’s can polynomial for each temporal be obtained by finding the roots of frequency followed by spectrum estimation for each source. The following result, which is an extension of the corresponding result in [8] to the broadband model under consideration, forms the basis of steps A) and B). be a block Toeplitz matrix of size Theorem 2: Let given as
, and the last equality where is some integer, Thus, it follows that follows from the fact that are orthogonal to those of Since for the columns of the full rank of , its columns span the columns orthogonal complement to the range space of the onto the subspace is equal to of , and the projection In the light of the above result, (13) can be rewritten as (23) where tr
(24)
tr
(25a)
tr
(25b)
As the trace is a linear operator, we get tr
(26)
Moreover, using the commutative property of the convolution operation, i.e., (27) where
is block data matrix defined as (28)
and each block is a circulant matrix
(28a) (20)
we get (29)
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
The Constraint Set As in [8], the ML solution for the vector is sought over an appropriate constraint set that is chosen to guarantee nontriviality as well as optimality. The nontriviality constraint (i.e., ensures that there is a unique correspondence between and for each and renders invertible since is block Toeplitz of whose each sub-block is a circulant matrix as given in (28). Typically, such a nontriviality constraint to be a monic polynomial is imposed by either requiring (henceforth referred to as the linear constraint) or in (henceforth called the norm constraint). setting It has been shown by Nagesha and Kay [15] in the narrowband context that estimates obtained using the norm constraint are better than those based on the linear constraint. Furthermore, since, in the array processing context, we are looking at the case where is a pure, complex sinusoid, it should have its roots on the unit circle in follows that This may be achieved by imposing conjugate the -plane polynomial , i.e., symmetry on the (30) This, of course, is only a necessary condition and not a sufficient one for the roots to lie on the unit circle. However, simulations suggest that this is quite adequate in most situations. In fact, for the case of a single source, this is known to be a necessary and sufficient condition. This fact has been used by Shaw [16] to guarantee all zeros to lie on the unit circle. This method can easily be adapted for use here. Finally, we recollect from our discussion in Section II that of each of the polynomials are related the roots Unfortunately, via (5), i.e., there is no simple way to express this requirement as an algebraic constraint on the coefficients of the polynomial It is because of this difficulty that we resort to the two-step procedure discussed in Section II and elaborated further in the next section to obtain a single (combined) from individually obtained estimates of estimate of each
623
to the wideband case, the IQML algorithm can be summarized as follows. ; choose an initial vector , and a) Initialization: Set denote it by b) Compute (32) is computed from in the same manner where in (20) from vector as c) Solve the quadratic minimization problem (33) over the appropriate constraint set (c.s.). by rooting , and d) Obtain the vector solve the minimization problem as given in (17) to This can be done by carrying out a simple, obtain 1-D search for each value of e) Construct the vector from f) Set . g) Check for convergence if yes go to step h) if not go to step b). . These are the desired spatial h) Find roots of ’s of the signals. parameters In practice, experience shows that there is not much to be lost by taking out step d) from the above iterative procedure and obtaining a least squares combined estimate of directly from the converged value of This fact can be used to advantage to reduce the computational complexity of the algorithm. Furthermore, simulations show that the least squares search itself can be replaced with a simpler operation of using for (34) where
V. A WIDEBAND ML ESTIMATION ALGORITHM
(34a)
The ML estimate of the 1-D predictor polynomials can be obtained by minimizing
without any discernible degradation in performance of the estimator. (31)
over an appropriate constraint set as discussed earlier. This constrained nonlinear minimization problem can be solved by adapting the IQML algorithm of [8], which requires the solution of a quadratic minimization problem at each step, and generally converges in a small number of steps.2 Adapted 2 A few words of caution are in order here when using the IQML. It has been pointed out that IQML may not converge to the deterministic maximum likelihood estimates [17]. In fact, IQML estimates are known to be biased and inconsistent. A modified IQML procedure has recently been suggested, in [18], that mitigates some of these deficiencies of the IQML procedure. It is indicated later in this section that it is easy to incorporate a similar modification in the broadband case to effect performance improvement over IQML.
Some Remarks on the Implementation of the Constraint Set As per the aforementioned discussion, we need to impose the constraints of conjugate symmetry and nontriviality on the It is possible to incorporate these conpolynomials straints in the iterative process by modifying the data matrix This can be easily done by making simple modification of the procedure suggested in [8]. We first describe the modified procedure for the conjugate symmetry constraint. (This is being done Let be an odd number s.t. is even, equivalent equations only for explanation. When s.t. can be developed.) Partition and (35a)
624
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
where
and (35b) For conjugate symmetry of that
is a matrix of size , given as of length .. .
, it can be easily seen (36)
, , and
is a vector
(44)
Combining both the constraints, the ML optimization problem reduces to a closed-form solution given by
where is a block permutation matrix of dimension and is given by
(45) Modified IQML (37)
In (37)
A modification of the IQML procedure has recently been suggested by Kristensson et al. [18], and it involves calculation of the signal subspace from the estimates obtained in the final has been modified to iteration. Cost function tr
(46)
(37a) where and is a matrix all zeros of size defined as
Let
and
be (38)
and (39) This reduces the minimization problem in (33) to
tr
(46a)
This modification has been shown to yield significant performance improvement in the narrowband case. While the details are omitted here, in the following, we also consider the impact of this modification for the wideband case. VI. SIMULATION RESULTS
(40) Define (41) It can then be shown [8] that solving the above minimization problem is equivalent to (42) Next, we consider the incorporation of the nontriviality constraints. As far as the norm constraint is constrained, it is easily implemented by searching for the eigenvector corresponding to the minimum eigenvalue of the matrix For implementation of the monicity constraint, we note that to be monic is equivalent to constraining constraining to equal unity As the last element of the polynomial our temporal frequencies are harmonically related, each coefis actually the DFT of the corresponding part ficient of Therefore, of vector at the frequency corresponding to the above constraint can be met easily by constraining the DFT of last elements of to be unity for all frequencies. In terms of the vector [which is defined by (38)], this constraint can be easily expressed as (43)
The approach proposed in this paper has been subjected to simulation experiments to study various issues such as 1) validity of the proposed method; 2) comparison of different formulations of the wideband IQML; 3) effect of deviations from the assumed harmonic model. For ease of demonstration, each source is described by a A ULA of sensors four frequency harmonic model with sensor separation equal to one half of the wavelength corresponding to the highest frequency is considered. All sources in each of the experiments are assumed to have the same power. A. Validity of the Proposed Method Simulation experiments have been carried out with a uniIn Fig. 1, we have form linear array of 24 sensors considered a scene with three uncorrelated sources in the field at and relative to the broadside of the array. For these relatively simple conditions, the norm and the linear constraints have almost similar performance. Furthermore, there appears to be little or no difference under these conditions between the IQML and the modified IQML. For this reason, results corresponding to the case of linear constraint only are presented here. The mean square error (MSE) in estimating the DOA’s for each of the sources for various SNR’s is computed by averaging the squared error over 100 independent trials. The results shown in Fig. 1 are for
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
(a)
625
(b)
(c) Fig. 1. Mean square error (MSE) and the CR bound “versus” SNR for three widely separated, wideband, uncorrelated sources (M = 24; N = 1): (a) Source at 20 . (b) Source at 45 . (c) Source at 80 :
(a)
(b)
(c) Fig. 2. MSE and the CRB “versus” SNR for three widely separated, wideband, fully correlated sources (M = 24; N = 1): (a) Source at 20 . (b) Source at 45 . (c) Source at 80 :
626
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
(a)
(b)
Fig. 3. Average MSE versus SNR for three wideband sources ( = 20 ; 45 ; 80 ), and M = 24. (a) For uncorrelated source case. (b) For correlated source case.
the case of a single snapshot and compared with the corresponding values of the Cram´er–Rao (CR) bound, which has been calculated by generalizing the results of [19]–[21] for the wideband formulation in this paper. As expected, the performance is seen to be quite close to the CR bound, especially for sources close to the broadside at high SNR. Fig. 2 shows the corresponding results when the three sources are fully correlated. The correlation of the sources is seen to have very little degrading effect on the performance except at low SNR’s. In each of the two cases (i.e. for correlated as well as uncorrelated cases), there appears to be a threshold SNR, below which the performance deviates significantly from the CR bound. Clearly, this effect is more pronounced for the case of correlated sources than otherwise. By increasing the number of snapshots, the number of unknown parameters to be estimated increases due to the increase in the number of signal sample values [7], [21]. However, the number of DOA parameters remains the same. It may be expected that the estimation of these direction parameters will improve with an increase in the number of snapshots. In order to study the behavior of estimators with a different number of snapshots, simulations have also been and for both the correlated carried out with and uncorrelated source cases. Fig. 3 summarizes this behavior by plotting the arithmetic average of the mean-square error of all the three sources. It is seen from this figure that the performance of the estimators indeed improves with increasing Finally, a case of two sources is considered with relative to the broad side of array, the first source fixed at whereas the bearing of the second source is varied. Fig. 4 compares the resulting MSE in the DOA estimates of the second source with the CR bound as its bearing is varied. and an These results are based on 50 snapshots SNR of 10 dB. As seen here, it appears to be possible to resolve even closely spaced sources with sufficient accuracy. B. Comparison of Different Formulations The preceding set of simulations have been conducted for relatively mild conditions when the number of sensors is large with only three well-separated sources in the field. As mentioned earlier, under these conditions, all the formulations
Fig. 4. MSE in DOA of a wideband source as a function of its separation from the other source present at 10 ; M = 24; N = 50; and SNR = 10 dB, with respective CRB.
Fig. 5. Comparison of IQML with linear constraint, IQML with the norm constraint, and MIQML. Average MSE of three wideband sources at 10 ; 20 ; and 30 versus SNR for M = 5; N = 10 000:
including MIQML yield similar results. In the next set of simulation studies, we make the conditions somewhat more stringent with a view to bringing out the differences, if any. First, consider the case of three sources present at and and a ULA of only five sensors Fig. 5 shows the arithmetic mean of the MSE in the DOA estimates of the three sources at various SNR’s. Simulations have carried These MSE’s are out with 10 000 snapshots average values obtained with 100 independent trials. It is seen that under these conditions, the IQML method with the norm constraint yields much better results than that with the linear constraint. Furthermore, at low SNR’s, the MIQML is seen to
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
627
ULA corresponding to the th snapshot can now be expressed as
(47)
Fig. 6. Comparison of IQML with linear constraint, IQML with the norm constraint, and MIQML. Average MSE of three wideband sources at 10 ; 20 ; and 30 versus snapshots for M = 5 and SNR = 10 dB.
where ’s are the damping constants corresponding to different frequency components. For nonzero value of , the model of (47) would constitute a significant deviation from the assumed harmonic model. The effect of this deviation is best brought out by studying the case of a single source. For convenience, the damping constant of each component is assumed to be the same. Fig. 8 shows the variation of the resulting MSE in the DOA estimates as the damping constant is varied. Fig. 8(a) shows the case of a 24-sensor array, 50 snapshots , and an SNR of 10 dB. Fig. 8(b) shows the results for a two sensor array 1000 snapshots and a 30-dB source. For both simulations, it is seen that all three methods exhibit a graceful degradation in performance as the value of the damping constant is increased. This shows that the approach proposed in this paper may be expected to perform well in practical applications where the harmonic model assumption may not held exactly. VII. CONCLUSION
Fig. 7. Comparison of IQML with linear constraint, IQML with the norm constraint, and MIQML. MSE for wideband source as a function of its separation from another source at 10 ; M = 5; N = 10 000; SNR = 10 dB.
outperform both of these. Fig. 6 similarly shows the arithmetic mean of the MSE in DOA estimates of the three sources as a function of number of snapshots. MIQML estimates are seen to improve consistently with the snapshot index. This may be attributed to the fact that the estimate of the covariance matrix and, hence, of the noise variance improves with increasing snapshots. Furthermore, the linear constraint formulation of the IQML seems to break down, whereas with the norm constraint, it still performs reliably at the SNR (10 dB) under consideration. for the Next, consider a scene of two sources : one being at and the other with same array a variable bearing. Fig. 7 compares the resulting MSE in the DOA estimates of second source with the CR bound, as its bearing is varied. Once again, the same behavior pattern is seen: IQML with the linear constraint failing completely and MIQML outperforming the IQML with the norm constraint, especially for closely spaced sources. C. Effect of Deviations from the Assumed Model Finally, we consider the important issue of robustness of the method proposed here when the underlying assumption of a harmonic signal model is violated. This is done here by considering the case where each source output consists of a damped sinusoids. Thus, the th time sample of sum of the signal from the th source impinging from direction (relative to the broad side of the array) at the th sensor of a
Using a harmonic model for broadband sources, a special ARMA model is formulated for the space time samples at an array of sensors from multiple broadband sources. This model is employed to formulate the estimation problem of DOA’s of broadband sources in terms of the coefficients of this special ARMA model. The problem so formulated is shown to admit adaption of IQML like procedures for its solution. The performance of the estimator is shown to be close to the CR bound. APPENDIX A OF THE EXPANDED FORM OF SPECIAL ARMA MODEL GIVEN BY (17) PROOF
Consider the 2-D data samples given by (A1) is the sum of 2-D sinusoid, i.e.,
where
(A2) where
and
are defined in (5). Define
by (A3)
First, we show that there exists a set of coefficients s.t. (A4)
628
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
(a)
(b)
Fig. 8. Effect of deviations from harmonicity. Mean square error versus “damping constant” of actual model. (a) Single source = 10 dB. (b) Single source M = 2; N = 10 000 SNR = 30 dB. All the three methods yield similar performance.
M
= 24; N = 50; SNR
(A9)
Defining
Substituting (A2) in (A4), we get
and
similar to
, i.e., (A11)
(A5) Equation (A1) can be rewritten as and Define equation can be rewritten as
(A12)
so that the above Substituting (A4) in (A12), we get
(A13)
(A6) Again, using (A12), we get
implying that (A7)
(A14)
satisfyThis proves that with coefficients ing (A7), (A4) will hold. We need to show, however, that (A7) 2-D linear has a nontrivial solution. (A7) corresponds to predictor polynomial equations, which can be neatly written in a vector form as
Hence, we get the expanded form of the special ARMA model as
(A8) where is the matrix defined as (A9), shown at the top of the given as page. is a vector of length (A10) and is a vector of all ones of length The matrix is seen to have a block Vandermonde kind of structure. Its rows are block Vandermonde and are distinct. Hence, it is nonsingular, and there exists a unique vector s.t. Equation (A4) is satisfied. Alternatively, it follows from (A7) that given the vector , the unknown frequencies can be obtained as the roots of the polynomial given in (19c).
(A15) where
equals unity. REFERENCES
[1] M. Wax, T. J. Shan, and T. Kailath. “ Spatio-temporal spectral analysis by eigenstructure methods,” IEEE Trans. Acoust., Speech, Signal Processing, vol ASSP-32, pp. 817–827, Aug 1984. [2] H. Wang and M. Kaveh, “Coherent signal subspace processing for detection and estimation of angle of arrival of multiple wideband sources,” IEEE Trans. Acoust., Speech, Signal Processing, vol ASSP-33, pp. 823–831, Aug. 1985. [3] K. M. Buckley and L. J. Griffiths, “Eigenstructure based broadband source location estimation,” in Proc. IEEE ICASSP, Tokyo, Japan, 1986, pp. 1869–1872. [4] G. Su and M. Morf, “Modal decomposition signal subspace algorithms,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-34, pp 585–602, June 1986.
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
[5] G. Bienvenu, “Eigensystem properties of the sampled space correlation matrix,” Proc. IEEE, ICASSP, Boston MA, 1983, pp. 332–335. [6] Y. Grenier, “Wideband source location through frequency-dependent modeling,” IEEE Trans. Signal Processing, vol 42, pp. 1087–1096, May 1994. [7] M. A. Doron, A. J. Weiss, and H. Messer, “Maximum likelihood direction finding of wideband sources,” IEEE Trans. Signal Processing, vol. 41, pp. 411–414, Jan 1993. [8] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter estimation of superimposed exponential signals in noise,” IEEE Trans. Acoust., Speech, Signal Processing, vol ASSP-34, pp 1081–1089, Oct. 1986. [9] M. P. Clark and L. L. Scharf, “Two dimensional modal analysis based on maximum likelihood,” IEEE Trans. Signal Processing, vol. 42, pp. 1443–1451, June 1994. [10] S. M. Kay, V. Nagesh, and J. Slisbury, “Broadband detection based on two-dimensional mixed autoregressive model,” IEEE Trans. Signal Processing, vol. 41, pp. 2413–2428, July 1993. [11] Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” Proc. IEEE ICASSP, 1991, pp. 3073–3076. [12] R. A. Monzingo and T. W. Miller, Introduction of Adaptive Arrays. New York: Wiley, 1980. [13] S. Haykin, Advances in Spectrum Analysis and Array Processing, vol. III. Englewood Cliffs, NJ: Prentice-Hall, 1995. [14] T. J. Ulrych and R. W. Clayton, “Time series modeling and maximum entropy,” Phys. Earth Planetary Interiors, vol. 12, pp. 188–200, Aug. 1976. [15] V. Nagesha and S. Kay, “On frequency estimation with the IQML algorithm,” IEEE Trans. Signal Processing, vol. 42, pp. 2509–2513, Sept. 1994. [16] A. K. Shaw, “Maximum likelihood estimation of multiple frequencies with constraints to guarantee unit circle roots,” IEEE Trans. Signal Processing, vol. 43, pp. 796–799, Mar. 1995. [17] P. Stoica, J. Li, and T. Soderstrom, “On the inconsistency of the IQML,” Signal Process., vol. 56, no. 2, pp. 185–190, Jan 1997. [18] M. Kristensson, M. Jansson, and B. Ottersten, “Modified IQML and a statistically efficient method for direction estimation without eigendecomposition,” in Proc. IEEE ICASSP, 1998. [19] S. F. Yau and Y. Bresler, “A compact Cramer-Rao bound expression for parametric estimation of superimposed signals,” IEEE Trans. Signal Processing, vol. 40, pp. 1226–1229, May 1992. [20] M. P. Clark, “Cramer Rao bound for two dimensional deterministic modal analysis,” in Proc. 27th Asilomar Conf. Signals Syst. Comput., Pacific Grove, CA, Nov. 1993, pp. 1151–1156. [21] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood and Cramer Rao bound,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 37, pp. 720–741, May 1989. [22] N. Chotikakamthorn, “A pre-filtering maximum likelihood approach to multiple source direction estimation,” Ph.D. dissertation, Dept. Elect. Eng., Imperial Coll. Sci. Technol. Med., Univ. London, London, U.K.
629
Monika Agrawal received the B.Tech. degree in electrical engineering and the M.Tech. degree in electronics and communications from the Regional Engineering College, Kurukshetra, India, in 1993 and 1995, respectively. She is currently a research scholar at the Indian Institute of Technology, Delhi, where she is pursuing the Ph.D. degree. Her research interests include digital signal processing and array processing and their applications to communications.
Surendra Prasad (SM’94) received the B.Tech. degree in electronics and electrical communication engineering from the Indian Institute of Technology, Kharagpur, in 1969 and the M.Tech. and Ph.D. degrees in electrical communication from the Indian Institute of Technology, New Delhi, in 1971 and 1974, respectively. He has been with the Indian Institute of Technology, New Delhi, since 1971, where he is presently a Professor of Electrical Engineering and Coordinator of the Research and Training Programme in Telematics. He was a Visiting Research Fellow at the Loughborough University of Technology, Loughborough, U.K., from 1976 to 1977, where he was involved in developing algorithms for adaptive array processing for HF arrays. He was also a Visiting Faculty Member at the Pennsylvania State University, University Park, from 1985 to 1986. His teaching and research interests include communications as well as statistical and digital signal processing. He has been a consultant to a number of government agencies as well as industry in these and related areas. Currently, he is engaged in research in various aspects of statistical signal processing, including, among others, single and multichannel ARMA filtering and modeling, high-resolution spectrum estimation and array processing, higher order spectrum estimation and their applications to problems in digital communication, underwater acoustics, and seismic data processing. He has published more than 70 papers in these areas in reputed journals. He has also edited a Special Issue of the Journal of I.E.T.E. in March–April 1989 in the area of statistical signal processing and a book entitled Signal Processing (New Delhi: I.E.T.E.). Dr. Prasad was the recipient of the Vikram Sarabhai Research Award in electronics and telecommunications in 1987, the Shanti Swarup Bhatnagar Award for engineering sciences in 1988, and the Om Prakash Bhashin Prize for Research in Electronics and Communications for 1994. He was a CoChairperson for the “Indo-U.S. Workshop in One and Two Dimensions” held in New Delhi in November 1989. He is a Fellow of the Indian National Academy of Engineering, the Indian National Science Academy, and the Indian Academy of Sciences.
619
DOA Estimation of Wideband Sources Using a Harmonic Source Model and Uniform Linear Array Monika Agrawal and Surendra Prasad, Senior Member, IEEE Abstract— We consider here the problem of estimation of DOA’s of multiple wideband sources incident on a uniform linear array (ULA) in the presence of spatially and temporally white Gaussian noise (WGN). The approach presented here builds up on the IQML algorithm suggested by Bresler and Macovski for the case of narrowband DOA estimation. It is shown that the concept of ARMA model for the observed data vector for the narrowband case can be generalized to model an appropriately stacked, space-time data vector obtained by combining the spacetime samples. The coefficients of the corresponding 2-D predictor polynomial can be used to represent the null subspace of the wideband array steering matrix, and rooting of the polynomial at each frequency, separately, gives the DOA estimates. These separate estimates at multiple frequencies are combined into a single DOA estimate in a least squares sense. This leads to the formulation of an IQML-like procedure for the spatial parameter estimation of wideband sources. Like its narrowband counterpart, the proposed approach is applicable to both noncoherent and coherent sources. The performance of the proposed method is studied via extensive computer simulations and by comparison with the CR bounds. Index Terms— Array signal processing, broadband sources, direction-of-arrival estimation, maximum likelihood techniques, special ARMA model, uniform linear array.
I. INTRODUCTION
T
HE PROBLEM of estimating the source parameters by processing data from a sensor array has received much attention in recent years. In so far as the estimation of the parameters of the narrowband sources is concerned, the theory is well established, and a large body of literature exists. The subspace-based methods have evolved into a powerful class of methods for direction-of-arrival (DOA) estimation of narrowband sources. These methods have also helped the formulation and solution of the maximum likelihood (ML) approach to the problem. For wideband sources, however, the literature is somewhat scanty, although a number of researchers have focussed their attention to this problem in recent years. Wax et al. [1], Wang and Kaveh [2], Buckley and Griffiths [3], and Morf and Su [4] are among the earlier researchers in this field who suggested some useful and interesting solutions for the problem of DOA estimation of wideband sources. Both incoherent aggregation suggested in [1] and coherent signal subspace (CSS) technique proposed in [2] represent in Manuscript received November 6, 1997; revised August 21, 1998. The associate editor coordinating the review of this paper and approving it for publication was Dr. Eric Moulines. The authors are with the Department of Electrical Engineering, Indian Institute of Technology, Delhi, New Delhi, India. Publisher Item Identifier S 1053-587X(99)01328-8.
some way adaption of narrowband processing techniques to the wideband case. Bienvenu [5] and Buckley and Griffiths [3] have shown how to divide the spatio-temporal observation space into signal and noise subspaces and have suggested search functions to estimate the DOA’s. Grenier [6] has generalized this notion to a so-called filter-spatial domain. Morf and Su [4] have suggested a modal decomposition of the signal subspace for estimating the DOA’s. In this paper, we are concerned with the solution of the wideband DOA estimation problem based on maximum likelihood principles. The problem of maximum likelihood estimation has, of course, been addressed earlier by Doron et al. [6], who have proposed several search functions to obtain the ML estimate of DOA for wideband sources. However, such an exhaustive search is computationally expensive and impractical for real-time applications. In this respect, the approach proposed by Bresler and Macovski [8] for the narrowband case, and a uniform linear array (ULA), is of particular interest, as it formulates an iterative method called the iterative quadratic maximum likelihood (IQML) to obtain the ML solution. The IQML method involves the creation of a Toeplitz matrix that spans the subspace orthogonal to the signal subspace and can be parameterized by the coefficients of a prediction polynomial, whose roots yield the DOA estimates. In effect, therefore, use of this parameterized model converts the problem of estimating the DOA’s to that of estimating the coefficients of this model from which DOA’s can be easily calculated. Clark and Scharf [9] have extended this work to twodimensional (2-D) damped harmonic signals in additive Gaussian noise which, in principle, can take care of the wideband case. However, as pointed out by them, it is not possible to directly generalize the one-dimensional (1-D) results to the 2-D case due to the lack of a fundamental theorem of algebra for 2-D polynomials. They suggested a solution to the problem by using a prediction polynomial in one of the two dimensions and a polynomial that interpolates between the two. This paper builds on the ideas available for the parameter estimation of multiple narrowband sources in WGN to develop pseudo-maximum likelihood parameter estimators for the wideband case by adapting the harmonic model of Kay et al. [10] for wideband sources. It is shown that the concept of an ARMA model for the observed data vector used in the narrowband case [8] can be generalized to model an appropriately stacked, space-time data vector obtained by combining the spatio-temporal samples. As in the narrowband case, it is shown that the coefficients of the resulting ARMA model can be used to represent the null subspace of the wideband
1053–587X/99$10.00 1999 IEEE
620
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
array steering matrix associated with the model of the stacked observations for the harmonic signal model considered here. The resulting 2-D predictor polynomial is rooted for each frequency separately to give the DOA estimates. This leads us to the formulation of an IQML-like procedure for the parameter estimation of wideband sources that is similar to that for the narrowband case suggested in [8]. Thus, the difficulties associated with the 2-D model analysis due to lack of a theorem of algebra for 2-D polynomials or with pairing the roots when the 2-D problem is addressed by breaking it into two 1-D problems (as in [11]), are avoided in this approach.1 The solution, however, falls short of providing a single ML estimate for each source. This is because the IQML search is conducted over an expanded signal space, which in turn leads to generation of separate DOA estimates for different frequencies. For applications where it would be more desirable to obtain a single DOA estimate for each wideband source, we suggest a modification of the IQML search whereby the separate estimates at multiple frequencies can be combined into a single estimate of the DOA in an appropriately defined least squares sense. It is seen that the procedure leads to very good estimates, as evidenced by comparison with the Cram´er–Rao bound. This paper has been organized as follows. In Section II, we formulate the problem and develop its mathematical model. In Section III, we develop the ML criterion in a parametric framework. Sections IV and V constitute the main body of our work, where we unfold our solution. Results and simulation studies are presented in Section VI. II. MODELING
AND
Let be the total number of wideband sources present. In this paper, we make the usual assumption of to be known. may be estimated using the methods given Alternatively, in [2]. Thus, at the th sensor, the th time sample of the output signal corresponding to th snapshot is given by (2) is spatially and temporally white Gaussian where such that noise of variance (3) Define
(4) where
is the number of sensors present in the ULA. Let (5)
be the temporal and spatial frequency variables, respectively. Using this notation and combining (1)–(3), we can write
(6) The above equation can be rewritten in a vector form as
PROBLEM FORMULATION
As mentioned, we model a broadband signal here as the sum of a finite number of monochromatic propagating plane waves. Each source is assumed to have a spectral support of Hence, the th time sample of the signal from the th (relative to the broad side source impinging from direction of the array) at the th sensor of a ULA corresponding to the th snapshot can be expressed as [10]
(7) where
(8) and
(1) where ’s are harmonically related frequencies representing is the the spectral components of the sources. In addition, Fourier coefficient of th signal corresponding to , the th , where frequency, array element spacing; velocity of propagation of plane waves; number of monochromatic waves used to model the broadband source. 1 It has been brought to our attention by one of the reviewers that there is some similarity of our approach to that of N. Chotikakamthorn, as reported in his doctoral work at the Department of Electrical Engineering, Imperial College of Science, Technology and Medicine, University of London [22]. However, at the time of submitting this paper, we did not have access to this dissertation; we were, therefore, unable to comment on the precise relationship of our work with that of [22].
(9) to be a block It is convenient to regard the vector may Vandermonde vector. A block Vandermonde matrix now be defined as (9a) , and are where the temporal and spatial parameters of the signal. Defining , the system equation (7) can be rewritten in a concise form as (10) ’s and In this formulation, the unknown parameters are ’s. Hence, our estimation problem should be formulated as follows.
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
Given snapshots of the data vector , estimate ’s and ’s, where ’s are embedded in the steering matrix , and ’s are the components of the vector However, this estimation would require the solution of a highly nonlinear parameters optimization problem in terms of the Alternatively, we can consider the problem of esticontains the mating the spatial vector , whose component required DOA information. Even here, ideally, the estimation problem should be posed in such a way that the components of the estimated are related to each other via (5), i.e., In the developments that follow in Sections III and IV, it is seen that there is no obvious algebraic procedure for imbedding such a constraint on the parameter vector It is convenient, therefore, to split this problem into two steps. In the first step, we consider the modified problem of independently at all relevant frequencies, i.e., estimating Using (5), this step would lead to for separate estimates of : one corresponding to each frequency in the harmonic model. In the second step, we exploit the implied structural relationship between by (5) to obtain a single estimate for DOA of each wideband source in a least squares sense. These two steps are considered in more detail in the sequel. III. MAXIMUM LIKELIHOOD CRITERION As discussed in the proceeding paragraph, the first step of the problem is concerned with estimation of the parameter vector More precisely, given snapshots of the data vector , estimate ’s and ’s, where ’s are embedded in the ’s are the components of the vector steering matrix , and The solution of this problem is taken up in this section. It has been shown that under WGN, ML estimators and least square estimators are equivalent [12]. Hence, ML estimate of the signal parameters can be obtained by solving the nonlinear least squares problem (11) is the Euclidean norm. This expression may be where further simplified by first substituting the least square estimate in terms of as of (12) The least square estimate of the can then be shown to be given by [13]
where vector
arg
arg
tr
(13)
arg
tr
(14)
or, equivalently arg
and are the projection matrices for projection Here, onto the column space of and onto its orthogonal complement, respectively, and are given by
(15)
where tr is the trace operator, and given by covariance matrix of
621
is the estimated
(16) being the number of snapshots of the observation vector Having obtained the extended spatial parameter vector (of ) by solving (13) or (14), our next step is to look dimension for a -dimensional vector of unknown directions that best fits the estimated vector in the least squares sense. In order to do this, we recall that for the th source, the components of the vector are given by Assuming that the ML estimation of would bring about a decoupling of the estimates of its subvectors from each other, a reasonable (although, in principle, not exact) formulation of the least squares problem to estimate can be stated as follows. For the estimated such that vector , find an angle (17) is minimum. Such a least squares fit helps to combine the independently obtained DOA estimates at the various frequencies into a single effective estimate for each source. In the next section, we consider first the issue of estimating the spatial vector via a characteristic polynomial represenThe tation of the broadband, spatio-temporal observations detailed algorithm that combines the two steps of estimating and then in an iterative search framework is taken up in Section V. IV. ML CRITERION AND THE MODEL PARAMETERS It has been shown in the literature for the narrowband case [8] that the output vector observed on a ULA in the presence of WGN obeys a special ARMA model and that the ML estimates of its parameters are directly related to those of its coefficients. In effect, use of this model converts the problem of estimating the DOA’s to that of estimating the special ARMA parameters from which the DOA’s can be easily derived. Taking a clue from this equivalence, we consider the possibility of formulating a similar, equivalent for the ARMA model for the spatio-temporal samples broadband case, as given by (6). In view of the spatio-temporal , it is proposed that we use an nature of the samples expanded form of the special ARMA model for these samples. More specifically, we establish the following results for the space-time vector Theorem 1: The space-time vector admits the ARMA representation given by (18) where the set of coefficients ’s are assumed to form a 2-D linear predictor polynomial given by (18a)
622
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
with and denoting the temporal and spatial complex frequency variables, respectively. Proof: Appendix A outlines a proof of this equivalent by adopting the methodology of Ulrych and description of Clayton [14]. The polynomial may be regarded as the linear predictor polynomial for the noiseless signal component of (2). This model can be considered to be a generalization of the Ulrych and Clayton model [14] for the case of 2-D harmonic of series. For convenience, we define the vectors and , respectively, as lengths
where each block is a circulant matrix of size
given as
(20a)
Then,
can be written as (21) be a column of and
Proof: Let to see that for
Then, it is easy
(19a) and
(22)
(19b) Here, elements of the vector polynomial for each as
form a 1-D linear predictor
(19c) is the th element of the vector It follows where ’s of this polynomial are known, that if the coefficients ’s can be estimated by rooting the spatial frequencies This result forms the basis of our proposed approach for the parameter estimation problem for wideband (more precisely multichromatic) plane waves in spatially and temporally white Gaussian noise. The method is similar to that for the narrowband case and consists of following steps. A) Establish a relation to represent the null space of the array steering matrix associated with the model of the space-time stacked vector , as defined by (10), in terms of the coefficients of the special ARMA model. B) Formulate the minimization problem in terms of a constrained quadratic minimization procedure for obtaining the ARMA coefficients Like the narrowband case, once the coefficients of the are estimated, the unknown DOA’s can polynomial for each temporal be obtained by finding the roots of frequency followed by spectrum estimation for each source. The following result, which is an extension of the corresponding result in [8] to the broadband model under consideration, forms the basis of steps A) and B). be a block Toeplitz matrix of size Theorem 2: Let given as
, and the last equality where is some integer, Thus, it follows that follows from the fact that are orthogonal to those of Since for the columns of the full rank of , its columns span the columns orthogonal complement to the range space of the onto the subspace is equal to of , and the projection In the light of the above result, (13) can be rewritten as (23) where tr
(24)
tr
(25a)
tr
(25b)
As the trace is a linear operator, we get tr
(26)
Moreover, using the commutative property of the convolution operation, i.e., (27) where
is block data matrix defined as (28)
and each block is a circulant matrix
(28a) (20)
we get (29)
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
The Constraint Set As in [8], the ML solution for the vector is sought over an appropriate constraint set that is chosen to guarantee nontriviality as well as optimality. The nontriviality constraint (i.e., ensures that there is a unique correspondence between and for each and renders invertible since is block Toeplitz of whose each sub-block is a circulant matrix as given in (28). Typically, such a nontriviality constraint to be a monic polynomial is imposed by either requiring (henceforth referred to as the linear constraint) or in (henceforth called the norm constraint). setting It has been shown by Nagesha and Kay [15] in the narrowband context that estimates obtained using the norm constraint are better than those based on the linear constraint. Furthermore, since, in the array processing context, we are looking at the case where is a pure, complex sinusoid, it should have its roots on the unit circle in follows that This may be achieved by imposing conjugate the -plane polynomial , i.e., symmetry on the (30) This, of course, is only a necessary condition and not a sufficient one for the roots to lie on the unit circle. However, simulations suggest that this is quite adequate in most situations. In fact, for the case of a single source, this is known to be a necessary and sufficient condition. This fact has been used by Shaw [16] to guarantee all zeros to lie on the unit circle. This method can easily be adapted for use here. Finally, we recollect from our discussion in Section II that of each of the polynomials are related the roots Unfortunately, via (5), i.e., there is no simple way to express this requirement as an algebraic constraint on the coefficients of the polynomial It is because of this difficulty that we resort to the two-step procedure discussed in Section II and elaborated further in the next section to obtain a single (combined) from individually obtained estimates of estimate of each
623
to the wideband case, the IQML algorithm can be summarized as follows. ; choose an initial vector , and a) Initialization: Set denote it by b) Compute (32) is computed from in the same manner where in (20) from vector as c) Solve the quadratic minimization problem (33) over the appropriate constraint set (c.s.). by rooting , and d) Obtain the vector solve the minimization problem as given in (17) to This can be done by carrying out a simple, obtain 1-D search for each value of e) Construct the vector from f) Set . g) Check for convergence if yes go to step h) if not go to step b). . These are the desired spatial h) Find roots of ’s of the signals. parameters In practice, experience shows that there is not much to be lost by taking out step d) from the above iterative procedure and obtaining a least squares combined estimate of directly from the converged value of This fact can be used to advantage to reduce the computational complexity of the algorithm. Furthermore, simulations show that the least squares search itself can be replaced with a simpler operation of using for (34) where
V. A WIDEBAND ML ESTIMATION ALGORITHM
(34a)
The ML estimate of the 1-D predictor polynomials can be obtained by minimizing
without any discernible degradation in performance of the estimator. (31)
over an appropriate constraint set as discussed earlier. This constrained nonlinear minimization problem can be solved by adapting the IQML algorithm of [8], which requires the solution of a quadratic minimization problem at each step, and generally converges in a small number of steps.2 Adapted 2 A few words of caution are in order here when using the IQML. It has been pointed out that IQML may not converge to the deterministic maximum likelihood estimates [17]. In fact, IQML estimates are known to be biased and inconsistent. A modified IQML procedure has recently been suggested, in [18], that mitigates some of these deficiencies of the IQML procedure. It is indicated later in this section that it is easy to incorporate a similar modification in the broadband case to effect performance improvement over IQML.
Some Remarks on the Implementation of the Constraint Set As per the aforementioned discussion, we need to impose the constraints of conjugate symmetry and nontriviality on the It is possible to incorporate these conpolynomials straints in the iterative process by modifying the data matrix This can be easily done by making simple modification of the procedure suggested in [8]. We first describe the modified procedure for the conjugate symmetry constraint. (This is being done Let be an odd number s.t. is even, equivalent equations only for explanation. When s.t. can be developed.) Partition and (35a)
624
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
where
and (35b) For conjugate symmetry of that
is a matrix of size , given as of length .. .
, it can be easily seen (36)
, , and
is a vector
(44)
Combining both the constraints, the ML optimization problem reduces to a closed-form solution given by
where is a block permutation matrix of dimension and is given by
(45) Modified IQML (37)
In (37)
A modification of the IQML procedure has recently been suggested by Kristensson et al. [18], and it involves calculation of the signal subspace from the estimates obtained in the final has been modified to iteration. Cost function tr
(46)
(37a) where and is a matrix all zeros of size defined as
Let
and
be (38)
and (39) This reduces the minimization problem in (33) to
tr
(46a)
This modification has been shown to yield significant performance improvement in the narrowband case. While the details are omitted here, in the following, we also consider the impact of this modification for the wideband case. VI. SIMULATION RESULTS
(40) Define (41) It can then be shown [8] that solving the above minimization problem is equivalent to (42) Next, we consider the incorporation of the nontriviality constraints. As far as the norm constraint is constrained, it is easily implemented by searching for the eigenvector corresponding to the minimum eigenvalue of the matrix For implementation of the monicity constraint, we note that to be monic is equivalent to constraining constraining to equal unity As the last element of the polynomial our temporal frequencies are harmonically related, each coefis actually the DFT of the corresponding part ficient of Therefore, of vector at the frequency corresponding to the above constraint can be met easily by constraining the DFT of last elements of to be unity for all frequencies. In terms of the vector [which is defined by (38)], this constraint can be easily expressed as (43)
The approach proposed in this paper has been subjected to simulation experiments to study various issues such as 1) validity of the proposed method; 2) comparison of different formulations of the wideband IQML; 3) effect of deviations from the assumed harmonic model. For ease of demonstration, each source is described by a A ULA of sensors four frequency harmonic model with sensor separation equal to one half of the wavelength corresponding to the highest frequency is considered. All sources in each of the experiments are assumed to have the same power. A. Validity of the Proposed Method Simulation experiments have been carried out with a uniIn Fig. 1, we have form linear array of 24 sensors considered a scene with three uncorrelated sources in the field at and relative to the broadside of the array. For these relatively simple conditions, the norm and the linear constraints have almost similar performance. Furthermore, there appears to be little or no difference under these conditions between the IQML and the modified IQML. For this reason, results corresponding to the case of linear constraint only are presented here. The mean square error (MSE) in estimating the DOA’s for each of the sources for various SNR’s is computed by averaging the squared error over 100 independent trials. The results shown in Fig. 1 are for
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
(a)
625
(b)
(c) Fig. 1. Mean square error (MSE) and the CR bound “versus” SNR for three widely separated, wideband, uncorrelated sources (M = 24; N = 1): (a) Source at 20 . (b) Source at 45 . (c) Source at 80 :
(a)
(b)
(c) Fig. 2. MSE and the CRB “versus” SNR for three widely separated, wideband, fully correlated sources (M = 24; N = 1): (a) Source at 20 . (b) Source at 45 . (c) Source at 80 :
626
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
(a)
(b)
Fig. 3. Average MSE versus SNR for three wideband sources ( = 20 ; 45 ; 80 ), and M = 24. (a) For uncorrelated source case. (b) For correlated source case.
the case of a single snapshot and compared with the corresponding values of the Cram´er–Rao (CR) bound, which has been calculated by generalizing the results of [19]–[21] for the wideband formulation in this paper. As expected, the performance is seen to be quite close to the CR bound, especially for sources close to the broadside at high SNR. Fig. 2 shows the corresponding results when the three sources are fully correlated. The correlation of the sources is seen to have very little degrading effect on the performance except at low SNR’s. In each of the two cases (i.e. for correlated as well as uncorrelated cases), there appears to be a threshold SNR, below which the performance deviates significantly from the CR bound. Clearly, this effect is more pronounced for the case of correlated sources than otherwise. By increasing the number of snapshots, the number of unknown parameters to be estimated increases due to the increase in the number of signal sample values [7], [21]. However, the number of DOA parameters remains the same. It may be expected that the estimation of these direction parameters will improve with an increase in the number of snapshots. In order to study the behavior of estimators with a different number of snapshots, simulations have also been and for both the correlated carried out with and uncorrelated source cases. Fig. 3 summarizes this behavior by plotting the arithmetic average of the mean-square error of all the three sources. It is seen from this figure that the performance of the estimators indeed improves with increasing Finally, a case of two sources is considered with relative to the broad side of array, the first source fixed at whereas the bearing of the second source is varied. Fig. 4 compares the resulting MSE in the DOA estimates of the second source with the CR bound as its bearing is varied. and an These results are based on 50 snapshots SNR of 10 dB. As seen here, it appears to be possible to resolve even closely spaced sources with sufficient accuracy. B. Comparison of Different Formulations The preceding set of simulations have been conducted for relatively mild conditions when the number of sensors is large with only three well-separated sources in the field. As mentioned earlier, under these conditions, all the formulations
Fig. 4. MSE in DOA of a wideband source as a function of its separation from the other source present at 10 ; M = 24; N = 50; and SNR = 10 dB, with respective CRB.
Fig. 5. Comparison of IQML with linear constraint, IQML with the norm constraint, and MIQML. Average MSE of three wideband sources at 10 ; 20 ; and 30 versus SNR for M = 5; N = 10 000:
including MIQML yield similar results. In the next set of simulation studies, we make the conditions somewhat more stringent with a view to bringing out the differences, if any. First, consider the case of three sources present at and and a ULA of only five sensors Fig. 5 shows the arithmetic mean of the MSE in the DOA estimates of the three sources at various SNR’s. Simulations have carried These MSE’s are out with 10 000 snapshots average values obtained with 100 independent trials. It is seen that under these conditions, the IQML method with the norm constraint yields much better results than that with the linear constraint. Furthermore, at low SNR’s, the MIQML is seen to
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
627
ULA corresponding to the th snapshot can now be expressed as
(47)
Fig. 6. Comparison of IQML with linear constraint, IQML with the norm constraint, and MIQML. Average MSE of three wideband sources at 10 ; 20 ; and 30 versus snapshots for M = 5 and SNR = 10 dB.
where ’s are the damping constants corresponding to different frequency components. For nonzero value of , the model of (47) would constitute a significant deviation from the assumed harmonic model. The effect of this deviation is best brought out by studying the case of a single source. For convenience, the damping constant of each component is assumed to be the same. Fig. 8 shows the variation of the resulting MSE in the DOA estimates as the damping constant is varied. Fig. 8(a) shows the case of a 24-sensor array, 50 snapshots , and an SNR of 10 dB. Fig. 8(b) shows the results for a two sensor array 1000 snapshots and a 30-dB source. For both simulations, it is seen that all three methods exhibit a graceful degradation in performance as the value of the damping constant is increased. This shows that the approach proposed in this paper may be expected to perform well in practical applications where the harmonic model assumption may not held exactly. VII. CONCLUSION
Fig. 7. Comparison of IQML with linear constraint, IQML with the norm constraint, and MIQML. MSE for wideband source as a function of its separation from another source at 10 ; M = 5; N = 10 000; SNR = 10 dB.
outperform both of these. Fig. 6 similarly shows the arithmetic mean of the MSE in DOA estimates of the three sources as a function of number of snapshots. MIQML estimates are seen to improve consistently with the snapshot index. This may be attributed to the fact that the estimate of the covariance matrix and, hence, of the noise variance improves with increasing snapshots. Furthermore, the linear constraint formulation of the IQML seems to break down, whereas with the norm constraint, it still performs reliably at the SNR (10 dB) under consideration. for the Next, consider a scene of two sources : one being at and the other with same array a variable bearing. Fig. 7 compares the resulting MSE in the DOA estimates of second source with the CR bound, as its bearing is varied. Once again, the same behavior pattern is seen: IQML with the linear constraint failing completely and MIQML outperforming the IQML with the norm constraint, especially for closely spaced sources. C. Effect of Deviations from the Assumed Model Finally, we consider the important issue of robustness of the method proposed here when the underlying assumption of a harmonic signal model is violated. This is done here by considering the case where each source output consists of a damped sinusoids. Thus, the th time sample of sum of the signal from the th source impinging from direction (relative to the broad side of the array) at the th sensor of a
Using a harmonic model for broadband sources, a special ARMA model is formulated for the space time samples at an array of sensors from multiple broadband sources. This model is employed to formulate the estimation problem of DOA’s of broadband sources in terms of the coefficients of this special ARMA model. The problem so formulated is shown to admit adaption of IQML like procedures for its solution. The performance of the estimator is shown to be close to the CR bound. APPENDIX A OF THE EXPANDED FORM OF SPECIAL ARMA MODEL GIVEN BY (17) PROOF
Consider the 2-D data samples given by (A1) is the sum of 2-D sinusoid, i.e.,
where
(A2) where
and
are defined in (5). Define
by (A3)
First, we show that there exists a set of coefficients s.t. (A4)
628
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 47, NO. 3, MARCH 1999
(a)
(b)
Fig. 8. Effect of deviations from harmonicity. Mean square error versus “damping constant” of actual model. (a) Single source = 10 dB. (b) Single source M = 2; N = 10 000 SNR = 30 dB. All the three methods yield similar performance.
M
= 24; N = 50; SNR
(A9)
Defining
Substituting (A2) in (A4), we get
and
similar to
, i.e., (A11)
(A5) Equation (A1) can be rewritten as and Define equation can be rewritten as
(A12)
so that the above Substituting (A4) in (A12), we get
(A13)
(A6) Again, using (A12), we get
implying that (A7)
(A14)
satisfyThis proves that with coefficients ing (A7), (A4) will hold. We need to show, however, that (A7) 2-D linear has a nontrivial solution. (A7) corresponds to predictor polynomial equations, which can be neatly written in a vector form as
Hence, we get the expanded form of the special ARMA model as
(A8) where is the matrix defined as (A9), shown at the top of the given as page. is a vector of length (A10) and is a vector of all ones of length The matrix is seen to have a block Vandermonde kind of structure. Its rows are block Vandermonde and are distinct. Hence, it is nonsingular, and there exists a unique vector s.t. Equation (A4) is satisfied. Alternatively, it follows from (A7) that given the vector , the unknown frequencies can be obtained as the roots of the polynomial given in (19c).
(A15) where
equals unity. REFERENCES
[1] M. Wax, T. J. Shan, and T. Kailath. “ Spatio-temporal spectral analysis by eigenstructure methods,” IEEE Trans. Acoust., Speech, Signal Processing, vol ASSP-32, pp. 817–827, Aug 1984. [2] H. Wang and M. Kaveh, “Coherent signal subspace processing for detection and estimation of angle of arrival of multiple wideband sources,” IEEE Trans. Acoust., Speech, Signal Processing, vol ASSP-33, pp. 823–831, Aug. 1985. [3] K. M. Buckley and L. J. Griffiths, “Eigenstructure based broadband source location estimation,” in Proc. IEEE ICASSP, Tokyo, Japan, 1986, pp. 1869–1872. [4] G. Su and M. Morf, “Modal decomposition signal subspace algorithms,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-34, pp 585–602, June 1986.
AGRAWAL AND PRASAD: DOA ESTIMATION OF WIDEBAND SOURCES USING A HARMONIC SOURCE MODEL
[5] G. Bienvenu, “Eigensystem properties of the sampled space correlation matrix,” Proc. IEEE, ICASSP, Boston MA, 1983, pp. 332–335. [6] Y. Grenier, “Wideband source location through frequency-dependent modeling,” IEEE Trans. Signal Processing, vol 42, pp. 1087–1096, May 1994. [7] M. A. Doron, A. J. Weiss, and H. Messer, “Maximum likelihood direction finding of wideband sources,” IEEE Trans. Signal Processing, vol. 41, pp. 411–414, Jan 1993. [8] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter estimation of superimposed exponential signals in noise,” IEEE Trans. Acoust., Speech, Signal Processing, vol ASSP-34, pp 1081–1089, Oct. 1986. [9] M. P. Clark and L. L. Scharf, “Two dimensional modal analysis based on maximum likelihood,” IEEE Trans. Signal Processing, vol. 42, pp. 1443–1451, June 1994. [10] S. M. Kay, V. Nagesh, and J. Slisbury, “Broadband detection based on two-dimensional mixed autoregressive model,” IEEE Trans. Signal Processing, vol. 41, pp. 2413–2428, July 1993. [11] Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” Proc. IEEE ICASSP, 1991, pp. 3073–3076. [12] R. A. Monzingo and T. W. Miller, Introduction of Adaptive Arrays. New York: Wiley, 1980. [13] S. Haykin, Advances in Spectrum Analysis and Array Processing, vol. III. Englewood Cliffs, NJ: Prentice-Hall, 1995. [14] T. J. Ulrych and R. W. Clayton, “Time series modeling and maximum entropy,” Phys. Earth Planetary Interiors, vol. 12, pp. 188–200, Aug. 1976. [15] V. Nagesha and S. Kay, “On frequency estimation with the IQML algorithm,” IEEE Trans. Signal Processing, vol. 42, pp. 2509–2513, Sept. 1994. [16] A. K. Shaw, “Maximum likelihood estimation of multiple frequencies with constraints to guarantee unit circle roots,” IEEE Trans. Signal Processing, vol. 43, pp. 796–799, Mar. 1995. [17] P. Stoica, J. Li, and T. Soderstrom, “On the inconsistency of the IQML,” Signal Process., vol. 56, no. 2, pp. 185–190, Jan 1997. [18] M. Kristensson, M. Jansson, and B. Ottersten, “Modified IQML and a statistically efficient method for direction estimation without eigendecomposition,” in Proc. IEEE ICASSP, 1998. [19] S. F. Yau and Y. Bresler, “A compact Cramer-Rao bound expression for parametric estimation of superimposed signals,” IEEE Trans. Signal Processing, vol. 40, pp. 1226–1229, May 1992. [20] M. P. Clark, “Cramer Rao bound for two dimensional deterministic modal analysis,” in Proc. 27th Asilomar Conf. Signals Syst. Comput., Pacific Grove, CA, Nov. 1993, pp. 1151–1156. [21] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood and Cramer Rao bound,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 37, pp. 720–741, May 1989. [22] N. Chotikakamthorn, “A pre-filtering maximum likelihood approach to multiple source direction estimation,” Ph.D. dissertation, Dept. Elect. Eng., Imperial Coll. Sci. Technol. Med., Univ. London, London, U.K.
629
Monika Agrawal received the B.Tech. degree in electrical engineering and the M.Tech. degree in electronics and communications from the Regional Engineering College, Kurukshetra, India, in 1993 and 1995, respectively. She is currently a research scholar at the Indian Institute of Technology, Delhi, where she is pursuing the Ph.D. degree. Her research interests include digital signal processing and array processing and their applications to communications.
Surendra Prasad (SM’94) received the B.Tech. degree in electronics and electrical communication engineering from the Indian Institute of Technology, Kharagpur, in 1969 and the M.Tech. and Ph.D. degrees in electrical communication from the Indian Institute of Technology, New Delhi, in 1971 and 1974, respectively. He has been with the Indian Institute of Technology, New Delhi, since 1971, where he is presently a Professor of Electrical Engineering and Coordinator of the Research and Training Programme in Telematics. He was a Visiting Research Fellow at the Loughborough University of Technology, Loughborough, U.K., from 1976 to 1977, where he was involved in developing algorithms for adaptive array processing for HF arrays. He was also a Visiting Faculty Member at the Pennsylvania State University, University Park, from 1985 to 1986. His teaching and research interests include communications as well as statistical and digital signal processing. He has been a consultant to a number of government agencies as well as industry in these and related areas. Currently, he is engaged in research in various aspects of statistical signal processing, including, among others, single and multichannel ARMA filtering and modeling, high-resolution spectrum estimation and array processing, higher order spectrum estimation and their applications to problems in digital communication, underwater acoustics, and seismic data processing. He has published more than 70 papers in these areas in reputed journals. He has also edited a Special Issue of the Journal of I.E.T.E. in March–April 1989 in the area of statistical signal processing and a book entitled Signal Processing (New Delhi: I.E.T.E.). Dr. Prasad was the recipient of the Vikram Sarabhai Research Award in electronics and telecommunications in 1987, the Shanti Swarup Bhatnagar Award for engineering sciences in 1988, and the Om Prakash Bhashin Prize for Research in Electronics and Communications for 1994. He was a CoChairperson for the “Indo-U.S. Workshop in One and Two Dimensions” held in New Delhi in November 1989. He is a Fellow of the Indian National Academy of Engineering, the Indian National Science Academy, and the Indian Academy of Sciences.
