- Email: [email protected]
Broadband ML estimation under model order uncertainty
ARTICLE IN PRESS Signal Processing 90 (2010) 1350–1356
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Broadband ML estimation under model order uncertainty c ¨ Pei-Jung Chung a,, Mats Viberg b, Christoph F. Mecklenbrauker a
Institute for Digital Communications, The University of Edinburgh, UK Department of Signals and Systems, Chalmers University of Technology, Sweden c Institute of Communications and RF Engineering, Vienna University of Technology, Austria b
a r t i c l e i n f o
abstract
Article history: Received 8 June 2009 Received in revised form 5 November 2009 Accepted 6 November 2009 ¨ Dedicated to J.F. Bohme on the occasion of his 70th birthday Available online 26 November 2009
The number of signals hidden in data plays a crucial role in array processing. When this information is not available, conventional approaches apply information theoretic criteria or multiple hypothesis tests to simultaneously estimate model order and parameter. These methods are usually computationally intensive, since ML estimates are required for a hierarchy of nested models. In this contribution, we propose a computationally efficient solution to avoid this full search procedure and address issues unique to the broadband case. Our max-search approach computes ML estimates only for the maximally hypothesized number of signals, and selects relevant components through hypothesis testing. Furthermore, we introduce a criterion based on the rank of the steering matrix to reduce indistinguishable components caused by overparameterization. Numerical experiments show that despite model order uncertainty, the proposed method achieves comparable estimation and detection accuracy as standard methods, but at much lower computational expense. & 2009 Elsevier B.V. All rights reserved.
Keywords: Direction of arrival Broadband signals Maximum likelihood estimation Overparameterized models Unknown number of signals
1. Introduction The problem of estimating direction of arrival (DOA) plays an important role in array processing and a variety of applications such as radar, sonar, wireless communications and geophysics. The maximum likelihood (ML) approach is characterized by excellent statistical properties and robustness against small sample numbers, signal coherence and closely located sources. The standard ML approach assumes the number of signals, m, to be known and maximizes the likelihood function over an m-dimensional parameter set Mm . When m is unknown, traditional methods jointly estimate the number of signals and parameters of interest. The popular information theoretic criterion based approach [16,19] and the multiple hypothesis testing procedure [2,4,6,15]
Corresponding author.
E-mail addresses: [email protected] (P.-J. Chung), [email protected] (M. Viberg), [email protected] ¨ (C.F. Mecklenbrauker). 0165-1684/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.11.013
all belong to this category. More specifically, these methods compute ML estimates for a series of nested models and select the one that best fits the underlying criteria. Since ML estimates for all candidate models are required, the computational complexity becomes very high if the number of candidate models is large. Therefore, avoiding this full-search procedure can be a crucial step to improve computational efficiency. In this work, we pursue a max-search direction finding algorithm for unknown numbers of signals. The main idea is to consider the maximally hypothesized model and select relevant components associated with true parameters through simple hypothesis tests. In the narrow band case, the test can be conducted with an F-distributed statistic [8]. For broadband signals, the distribution under the null hypothesis has no closed form expression. Therefore, we apply the Cornish–Fisher expansion [10] to approximate the test threshold. In the previous publication [9], we exploited the important property that the overparameterized ML estimate contains the true parameters. This property was verified for narrow band, stochastic signal models [6]. The proof for the broadband,
ARTICLE IN PRESS P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
deterministic case will be presented here. Moreover, we shall discuss the identification problem caused by overparameterization and introduce a criterion to reduce indistinguishable components. The proposed procedure is different from the aforementioned model order estimation methods [6,14,17,19] in that it is aimed at extracting useful information about DOA parameters regardless of whether the number of signals is correctly specified or not. However, as a byproduct, the number of relevant components can be considered as an estimate for the actual number of signals. As it will be shown later, the proposed maxsearch procedure is computationally more attractive than the full search methods. The robustness against model order uncertainty could be of great importance to channel identification in wireless communications. In a time-varying propagation environment, our algorithm may help to distinguish the relevant point-like multipath components from the spurious ones with a minimum prior knowledge on model order. In the following, we give a brief description of the broadband signal model. Section 3 discusses asymptotic properties of ML estimation with overparameterized models. In Section 4, we develop the broadband ML estimation procedure and derive the Cornish–Fisher approximation to the test threshold. Simulation results are presented in Section 5, whereas Section 6 concludes the paper.
2. Broadband signal model Consider an array of n sensors receiving m broadband signals emitted by far-field sources located at h ¼ ½y1 ; . . . ; ym T . The sensor array outputs xðkÞ ðtÞ; ðt ¼ 0; . . . ; T1Þ are divided into K snapshots of length T 0 ¼ T=K and short-time Fourier-transformed T 1 1 X X ðkÞ ðoÞ ¼ pffiffiffi wðtÞxðkÞ ðtÞejot ; Tt¼0
ð1Þ
where fwðtÞgT1 t ¼ 0 is a window function. For large sample size, T 0 -1, the frequency domain data are approximately described by the following relation: X ðkÞ ðoÞ ¼ H m ðo; hm ÞS ðkÞ ðoÞ þ U ðkÞ ðoÞ;
ð2Þ
where the matrix H m ðo; hm Þ ¼ ½dðo; y1 Þ dðo; yi Þ dðo; ym Þ 2 Cnm consists of m steering vectors with the i th column dðo; yi Þ, corresponding to the i th incoming wave. We assume that the array manifold is unambiguous, i.e. any collection of m steering vectors corresponding to distinct parameters is linearly independent, and smooth, i.e. the third derivatives of dðo; yi Þ with respect to yi are bounded. The signal waveform S ðkÞ ðoÞ is considered as unknown and deterministic. The spatially white noise U ðkÞ ðoÞ results only from sensors. Under certain regularity conditions [3], X ðkÞ ðoj Þ ðk ¼ 1; . . . ; K; j ¼ 1; . . . ; JÞ are asymptotically independent, complex normally distributed with mean H m ðoj ; hm ÞS ðkÞ ðoj Þ and covariance matrix nðoj ÞI,
1351
where nðoj Þ is the unknown noise spectral parameter and I is an identity matrix of conformable dimension. Given a pre-specified number of signals, m, the broadband ML estimate is obtained by minimizing the negative log-likelihood function over the m-dimensional parameter space Mm [1,13]:
h^ m ¼ arg min lðhm Þ; hm 2Mm
lðhm Þ ¼
J X
logtr½ðIPðoj ; hm ÞÞC^ x ðoj Þ;
ð3Þ
j¼1
where Pðoj ; hm Þ represents the projection matrix onto the subspace spanned by the columns of H m ðoj ; hm Þ. The power spectral estimate of array outputs is give by P C^ x ðoj Þ ¼ ð1=KÞ Kk ¼ 1 X ðkÞ ðoj ÞX ðkÞ ðoj ÞH . In many array processing applications, the number of signals m is unknown. The most successful approaches exploit the information theoretic criteria or carry out multiple hypothesis testing to jointly estimate the DOA parameters and the model order. Given a hierarchy of nested models represented by nested parameter spaces Mm , i.e. hm 2 Mm for 1 r m rM, M 1 M2 MM ;
ð4Þ
these methods derive ML estimates for all candidate models and select the best one according to the underlying decision rule. Such full search procedures are typically computationally demanding, due to the multidimensional nonlinear minimization required by each candidate model. In contrary to this full search idea, we suggest a computationally more efficient alternative that considers the maximally hypothesized model MM and extracts DOA information from the corresponding estimate. 3. ML estimation of overparameterized models The key idea behind the proposed algorithm is that asymptotically, for an overparameterized model, the ML estimate contains components of the true DOA parameter, ho . By overparameterization, we mean that the assumed number of signals m is larger than the true number of signals mo . In other words, the asymptotic estimate derived from m 4mo ; hm ¼ ½y1 ; y2 ; . . . ; ym T , contains mo components associated with those of ho ¼ ½yo1 ; yo2 ; . . . ; yomo T . The remaining ðmmo Þ components correspond to redundant parameters. Applying the asymptotic theory of quasi-ML estimation [18], the analysis in [7] shows that under stochastic signal models, the narrow band ML estimator converges to the minimizing point of the Kullback–Leibler distance between the pdf of the assumed model and that of the true model. This overparameterized estimate hm is characterized by the following asymptotic property (valid for large K): spðH m ðo; hm ÞÞ*spðH mo ðo; ho ÞÞ;
ð5Þ
where spðH m ðo; hm ÞÞ and spðH mo ðo; ho ÞÞ denote the signal subspaces corresponding to hm and ho , respectively.
ARTICLE IN PRESS 1352
P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
Clearly, hm with the true parameter ho as a subvector is a solution to (5). These results can be generalized to broadband ML estimation with the deterministic (or conditional) signal model. However, rather than minimizing the Kullback– Leibler distance as in the stochastic case, we shall minimize the asymptotic likelihood function directly. This approach does not require standard assumptions on ML estimation such as fixed number of unknown parameters and i.i.d. samples which are not fulfilled in the current signal model. To investigate large sample properties, we shall exploit the true covariance matrix of array outputs H
C xo ðoÞ ¼ H mo ðo; ho ÞC so ðoÞH mo ðo; ho Þ þ no ðoÞI;
ð6Þ
where C so ðoÞ ¼ E½S ðkÞ ðoÞS ðkÞ ðoÞH and no ðoÞ represent the signal covariance matrix and noise level, respectively. It is well established that for K-1, the sample covariance matrix C^ x ðoÞ converges with probability one to the ensemble covariance matrix C xo ðoÞ. Result 1. Let m denote the assumed number of signals and m4 mo . For large sample size, the ML estimate obtained by minimizing the broadband likelihood function (3) converges to hm that minimizes the criterion Q ðhm Þ ¼
J X
logtr½ðIPðoj ; hm ÞÞC xo ðoj Þ
ð7Þ
j¼1
and achieves the minimum Q ðhm Þ ¼
J X
lognðoj Þ þ JlogðnmÞ:
ð8Þ
j¼1 j ; hm ÞÞ
Furthermore, the assumed signal subspace spðH m ðo and the true signal subspace spðH mo ðoj ; ho ÞÞ are asymptotically related as follows (valid for large K): spðH m ðoj ; hm ÞÞ*spðH mo ðoj ; ho ÞÞ;
j ¼ 1; . . . ; J:
ð9Þ
Proof. Recall that the finite sample likelihood function (3) is given by lðhm Þ ¼
J X
logtr½ðIPðoj ; hm ÞÞC^ x ðoj Þ:
rewritten as Q ðhm Þ ¼
J X
Q ðhm ; oj Þ;
j¼1
Q ðhm ; oj Þ ¼ logtr ½U 2 ðoj ÞH C xo ðoj ÞU 2 ðoj Þ:
ð10Þ
According to Corollary 4.3.18 of [11], the minimum of Q ðhm ; oj Þ subject to the constraint U 2 ðoj ÞH U 2 ðoj Þ ¼ I nm Q ðhm ; oj Þ ¼ min logtr ½U 2 ðoj ÞH C xo ðoj ÞU 2 ðoj Þ U 2 ðoj Þ
¼ log½ðnmÞnðoj Þ
ð11Þ
is achieved if the columns of the minimizing matrix U 2 ðoj Þ are orthonormal eigenvectors corresponding to the ðnmÞ smallest eigenvalues of C xo ðoj Þ. The ðnmÞ smallest eigenvalues are equal to the noise spectral parameter nðoj Þ because ðnmÞ oðnmo Þ. This also implies that spðU 2 ðoj ÞÞ is a subspace of the noise subspace of C xo ðoj Þ. Since U ðoj Þ ¼ ½U 1 ðoj Þ U 2 ðoj Þ forms an orthonormal basis for Cn1 , spðU 1 ðoj ÞÞ*spðH mo ðoj ; ho ÞÞ:
ð12Þ
The relation (9) is a direct consequence of spðU 1 ðoj ÞÞ ¼ spðHðoj ; hm ÞÞ and it is satisfied for j ¼ 1; . . . ; J if hm includes the true parameter ho as a subvector. As all frequency bins are minimized at hm , the broadband criterion Q ðhm Þ achieves its minimum (8) at hm . The minimum asymptotic likelihood is given by Q ðhm Þ ¼
J X
lognðoj Þ þ JlogðnmÞ:
j¼1
This completes the proof.
&
Eq. (8) shows the minimum asymptotic likelihood Q ðhm Þ varies with the assumed number of m. In contrast to the stochastic model whose asymptotic likelihood remains constant for m 4mo [7], the selection of relevant estimates can no longer rely on each component’s contribution to the likelihood function. Therefore, we propose a solution based on hypothesis testing in the following section.
j¼1
Following the fact C^ x ðoj Þ-C xo ðoj Þ for K-1, the asymptotic likelihood function becomes Eq. (7). Provided the array response (steering vector) is a smooth function of hm , it is straightforward to verify that the convergence is uniform in the parameters. Consequently, the parameter estimates converge to the minimizing arguments of (7) in the same statistical sense (e.g. w.p. 1) as the sample covariance matrix converges. To find the minimum of Q ðhm Þ, consider the singular value decomposition H m ðoj Þ ¼ Uðoj ÞRðoj ÞVðoj ÞH where Uðoj Þ ¼ ½U 1 ðoj Þ U 2 ðoj Þ is a unitary matrix. The n m matrix U 1 ðoj Þ consists of the first m columns of Uðoj Þ corresponding to the m largest singular values. Then we can express the projection matrix as Pðoj Þ ¼ U 1 ðoj Þ U 1 ðoj ÞH and its orthogonal complement as IPðoj Þ ¼ U 2 ðoj ÞU 2 ðoj ÞH . The asymptotic likelihood (7) can be
4. The max-search approach Motivated by the above observation, we suggest to compute the ML estimate for the maximally hypothesized model, m ¼ M, and select relevant components as estimates. More precisely, the proposed algorithm considers the estimate derived from an M-dimensional space
h^ M ¼ argminlðhM Þ; hM
ð13Þ
where M Zmo . The M 1 vector h^ M ¼ ½y^ 1 ; . . . ; y^ M T contains more elements than the mo 1 vector ho . The elements of h^ M that are associated with those of ho , are referred to as relevant components. The remaining ðMmo Þ components of h^ M are the redundant components. Clearly, the key step in the proposed max-search procedure is the identification of relevant components.
ARTICLE IN PRESS P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
In the stochastic model, relevant components are selected by thresholding the likelihood function [5], because redundant components do not contribute to the likelihood function. For the deterministic case, the asymptotic minimum (8) varies with assumed model order m and does not allow such simple selection scheme. We suggest the following hypothesis test to validate the ith component:
AM : XðoÞ ¼ HM ðo; h^ M ÞS M ðoÞ þUðoÞ;
ð15Þ
j¼1
J 1X n1 log 1 þ F i ð oj Þ : n2 J j¼1
ð16Þ
Here, Fi ðoj Þ can be interpreted as an incremental signalto-noise ratio F i ð oj Þ ¼
n2 tr½ðPðoj ; h^ M ÞPðoj ; h~ i ÞÞC^ x ðoj Þ n1 tr½ðIPðoj ; h^ M ÞÞC^ x ðoj Þ
and it turns out to be asymptotically Fn1 ;n2 distributed (valid for large K) under suitable assumptions (see Section 4.1). The associated degrees of freedom are [14] n1 ¼ 3K;
n2 ¼ Kð2n2M1Þ:
ð17Þ
The component y^ i is relevant if Hi is rejected. Given a significance level a, the null hypothesis Hi is rejected if Ti exceeds a threshold ta . Consequently,
y^ i is relevant if Ti Z ta :
ð18Þ
The output of the algorithm is the relevant vector that contains all components
h^ o ¼ ½y^ ð1Þ ; . . . ; y^ ðpÞ T :
(2) Compute test statistic Ti , i ¼ 1; . . . ; M. (3) Test Hi against AM , i ¼ 1; . . . ; M. y^ is relevant if T Z ta . i
ð14Þ
contains all elements of h^ M except the ith component. The n ðM1Þ matrix H M1 ðh~ i Þ contains steering vectors corresponding to the DOA parameters in h^ M1 . Given the estimate h^ M , the hypothesis test (14) decides whether the signal Si ðoÞ associated with the ith component is zero. Note that Si ðoÞ ¼ 0 implies y^ i does not correspond to any actual signal source. The test statistic Ti ði ¼ 1; 2; . . . ; MÞ derived from the likelihood ratio principle, lT ðh~ i ÞlT ðh^ M Þ, is given by ! J tr½ðIPðoj ; h~ i ÞÞC^ x ðoj Þ 1X Ti ¼ log J tr½ðIPðoj ; h^ M ÞÞC^ x ðoj Þ ¼
Input: fX ðkÞ ðoj Þ: k ¼ 1; . . . ; K; j ¼ 1; . . . ; Jg, upper bound on number of signals: M, test threshold: ta . (1) Find the ML estimate h^ M .
Output: relevant vector h^ o ¼ ½y^ ð1Þ ; y^ ð2Þ ; . . . ; y^ ðpÞ T .
where Hi and AM represent the null hypothesis and the alternative, respectively. The ðM1Þ 1 vector
h~ i ¼ ½y^ 1 y^ i1 y^ i þ 1 y^ M T
Table 1 The max-search algorithm for broadband ML estimation with unknown numbers of signals.
i
Hi : XðoÞ ¼ HM1 ðo; h~ i ÞS~ M1 ðoÞ þ UðoÞ;
1353
ð19Þ
As a byproduct of the algorithm, the number of relevant components provides an estimate for the number of signals. However, we emphasize that the primary concern of the proposed approach is parameter estimation. Unlike methods designed for model order determination, whose performance is measured solely by correctness of the number of signals, the output of the proposed algorithm h^ o still contains useful information about the true parameters ho even when the number of signals is misspecified. Given an upper bound on the
number of the signals M and the threshold ta , the proposed max-search algorithm is summarized in Table 1.
4.1. Cornish–Fisher approximation Under the null hypothesis Hi , the statistic Fi ðoj Þ is asymptotically Fn1 ;n2 distributed with the degrees of freedom n1 , n2 given by (17) for large K. As (20) shows, in the narrow band case, J ¼ 1, the test can be equivalently conducted with an F-statistic. However, for broadband signals, the distribution under the null hypothesis has no analytical expression. To overcome this difficulty, we suggest the Cornish–Fisher expansion to approximate the test threshold. The summands in (16) can be interpreted as i.i.d. samples drawn from the random variable n1 ð20Þ Fi ; T i ¼ log 1 þ n2 where the statistic F i is Fn1 ;n2 - distributed. This property follows from the asymptotic independence of the Fouriertransformed data X ðkÞ ðoj Þ ðj ¼ 1; . . . ; JÞ and the same degrees of freedom, n1 and n2 , for each frequency bin. Let T i ¼ J 1=2
ðT i mT Þ
sT
ð21Þ
denote the standardized statistic where mT ¼ E½T i and s2T ¼ E½T i mT 2 are the mean and variance of T i , respectively. With the knowledge of mean and variance, the distribution of T i can be approximated by the Edgeworth expansion [10]. For a statistic whose distribution admits an Edgeworth expansion, the a- quantile, defined by PfT i rt a g ¼ 1a
ð22Þ
can be approximated by the Cornish–Fisher expansion, 1 1 1 t a ¼ za þ pffiffi p1 ðza Þ þ p2 ðza Þ þ þ r=2 pr ðza Þ þ oðJr=2 Þ; J J J ð23Þ where za is the standard normal aquantile. The polynomial pr ðzÞ is of degree at most r þ 1, and depend on the cumulants only up to order r þ2. This expansion is valid for fixed r, as J-1. The quality of the approximation is given by the remainder oðJr=2 Þ. For fixed samples, the accuracy of the approximation can be improved by including more terms in the expansion. The first two
ARTICLE IN PRESS 1354
P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
polynomials are given by p1 ðzÞ ¼
k3 6
Sample mean of relevant components 37
ðz2 1Þ;
36
k 3 2 k 4 k4 2 z ðz 1Þ þ z ðz 3Þ þ ðz 10z þ 15Þ : p2 ðzÞ ¼ 18 24 72 2 3
35
r
where Cðr1Þ ðsÞ ¼ dr =ds ðlogGðsÞÞ denotes the r th derivative of the logarithm of the gamma function. The mean mT and variance s2T are related to the cumulants as follows:
s2T ¼ k2 :
33 32 31 30 29
ð25Þ
28
Computational details of cumulants kr are presented in Appendix.
27
mT ¼ k1 ;
θ(1), M = 4 θ(2), M = 4 θ(1), M = 3 θ(2), M = 3
34
The r th cumulant kr of T i depends on the degrees of freedom n1 , n2 through the following formula: h n n n i kr ¼ ð1Þr Cðr1Þ 1 Cðr1Þ 1 þ 2 ; ð24Þ 2 2 2
Mean (degree)
2 3
−10 −8
−6
−4
−2
0
2
4
6
8
10
SNR (dB)
4.2. Identification of indistinguishable components One situation that may occur in the proposed algorithm is that the estimate h^ M contains indistinguishable components, which leads to a rank deficient steering matrix Hðo; h^ M Þ. By ‘‘indistinguishable components’’ we mean that two DOA parameters are too close or almost identical at the minimum of the likelihood function. Under the assumption of unambiguous array manifold, property (9) permits such parameters because the underlying model is overparameterized. In this case, both the relevant vector h^ o and the test (14) no longer provide satisfying results. To overcome difficulties caused by overparameterization, we suggest to examine the smallest singular value s^ M of the steering matrix Hðo; h^ M Þ.1 If js^ M j r d where d is a pre-specified small positive number, then the procedure computes the ML estimate for the next reduced model MM1 . If the resulting Hðo; h^ M1 Þ is full rank, then the relevant components are chosen from h^ M1 . Otherwise, the same step is repeated for the next reduced model. In the simulation, we observe that this modified step is carried out only by 5% of all trials. The computational cost is still much lower than the full search approach. The indistinguishable components are usually almost identical; for example, the second and third components of the h^ 4 ¼ ½37:93783 29:21293 29:21353 60:39263 estimate differ from each other by less than 0:0063 . Hence, it is sufficient to consider the steering matrix at one frequency bin. In Section 5, the threshold d is chosen based on numerical data. How to select d in a statistically justified manner is still under investigation. 5. Simulation In the simulation, a uniform linear array of 10 sensors with inter-element spacings of half a wavelength l=2 is employed. The wavelength is given by l ¼ v=f0 where v represents the propagation velocity and f0 is a pre-selected 1 Alternatively, a minimum separation in the elements of h^ M could be used.
Fig. 1. Sample mean of relevant components. M ¼ 3, 4. The true DOA parameter ho ¼ ½283 363 , SNR ¼ ½10 : 2 : 10 dB.
reference frequency. The broadband signals are generated by m0 ¼ 2 uncorrelated signals, located at h0 ¼ ½283 363 with respect to array broadside, and of various strengths. We choose J ¼ 11 frequency bins equally spaced between ð17=32Þf0 and f0 for processing. The number of snapshots is K ¼ 50. The difference in signal strengths is [1 0] dB, where 0 dB corresponds to the reference signal. The signal to noise ratio (SNR), defined as 10 log10 E½jSi ðoj Þj2 =nðoj Þ for the i th signal ði ¼ 1; 2Þ, varies from 10 to 10 dB in a 2 dB step. We consider two upper bounds on the number of signals: M ¼ 3 and 4. The latter represents a larger mismatch in the model order. Each experiment consists of 200 Monte Carlo trials. The test level is chosen to be a ¼ 0:05. The parameter d ¼ 0:02 is used to test whether the steering matrix is rank deficient or not. Fig. 1 shows the sample mean of the relevant components y^ ð1Þ , y^ ð2Þ , respectively. For both M ¼ 3 and 4, the bias is less than 0.061 over the entire SNR range. Since the results are obtained from finite samples, we conjecture that y^ ð1Þ and y^ ð2Þ are asymptotically bias free. The empirical variance of y^ ð1Þ and y^ ð2Þ is presented in Fig. 2. For comparison, we also show results obtained from the no-mismatch case with M ¼ m0 ¼ 2. All three curves decline with increasing SNR. The estimates obtained from M ¼ m0 ¼ 2 have the smallest variance. For both relevant components y^ ð1Þ and y^ ð2Þ , M ¼ 4 results in a larger variance than M ¼ 3. This suggests that the variance increases with a larger degree of mismatch. In addition, it is observed in the simulation that bias and variance can be reduced by increasing the number of data samples K or frequency bins J. In Fig. 3, we compare the probability of correct detection of the proposed algorithm with the multiple hypothesis testing procedure [6]. By ‘‘correct detection’’, we mean that the number of relevant components equals the true number of signals. At the low SNR region, 10 to 4 dB, the curves below 100% indicate that the signals are not always detected. However, both M ¼ 3 and 4 have a higher probability of correct detection than the multiple test. For SNR 6 to 10 dB,
ARTICLE IN PRESS P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
Empirical variance of θ(1)
101 Variance (degree2)
1355
100
θ(1), M = 4
10−1
θ(1), M = 3 θ(1), m0 = 2
10−2 −10
−8
−6
−4
0 SNR (dB)
2
4
6
8
10
4
6
8
10
Empirical variance of θ(2)
101 Variance (degree2)
−2
100 θ(2), M = 4
10−1
θ(2), M = 3 θ(2), m0 = 2
10−2 −10
−8
−6
−4
−2
0
2
SNR (dB) Fig. 2. Empirical variance of relevant estimates y^ ðiÞ , i ¼ 1, 2.
signals is unknown or misspecified. The price paid for model uncertainty is increased variance. Furthermore, the comparison between the proposed algorithm and the full search procedure [6] showed that without computing ML estimates for each candidate model, the former has a lower SNR threshold and achieves good detection performance at high SNRs. The computational time required by the suggested max-search procedure is on average 35% less than the full search procedure. The saving is expected to grow for larger numbers of candidate models.
Estimation of Number of Signals 1 Probability of Correct Detection
0.9 0.8
M=4 M=3 Multiple test
0.7 0.6 0.5 0.4 0.3
6. Conclusion
0.2 0.1 0 −10 −8
−6
−4
−2
0
2
4
6
8
10
SNR (dB) Fig. 3. Probability of correct detection. Comparison of multiple test procedure [5], M ¼ 3, 4.
the multiple test procedure achieves 100% probability of correct detection, while M ¼ 3; 4 increase from 90% to 95%. The curve associated with M ¼ 3 shows a slightly higher probability of correct detection than that of M ¼ 4. In summary, simulation results demonstrate that relevant components convey useful information about the true parameters even when the correct number of
We developed a broadband ML estimation procedure for unknown numbers of signals. The suggested algorithm computes ML estimates only for the maximally hypothesized number of signals. The relevant components associated with the true DOA parameters are selected by simple hypothesis tests. The test threshold is approximated by the Cornish–Fisher expansion. We also introduce a criterion to reduce indistinguishable components caused by overparameterization. Compared to traditional methods for joint parameter estimation and signal detection, the proposed max-search approach avoids the full search process through a series of nested models, which leads to significant improvement in computational efficiency if the number of candidate models is large. Numerical results showed that the proposed algorithm achieves comparable estimation accuracy
ARTICLE IN PRESS 1356
P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
as the standard ML approach does. The number of signals can also be accurately determined by the number of relevant components. The proposed algorithm provides a computationally attractive alternative to existing methods, particularly in the low SNR region. A potential drawback of the proposed approach is that over-specification of the number of signals may lead to higher sensitivity to errors in the steering vector model. This is an interesting topic for further investigations.
The r th derivative of (30) leads to the desired r th cumulant kr ðr ¼ 1; 2; . . .Þ of T dr kr ¼ r logfT ðsÞ ds s¼0 r d dr ¼ log G ðqjsÞ log G ðqjs þ pÞ r r ds ds s¼0 ¼ ð1Þr ½Cðr1Þ ðqÞCðr1Þ ðq þ pÞ:
ð31Þ
We obtain (24) by replacing p ¼ n1 =2, q ¼ n2 =2 into the above equation.
Acknowledgements
References
The authors wish to thank the anonymous reviewers for their comments that significantly improved the manuscript and also thank the Managing Editor Prof. A.M. Zoubir for coordinating a speedy review.
¨ [1] J.F. Bohme, Array processing, in: S. Haykin (Ed.), Advances in Spectrum Analysis and Array Processing, Prentice-Hall, Englewood Cliffs, NJ, 1991, pp. 1–63. ¨ [2] J. F. Bohme, Statistical array processing of measured sonar and seismic data, in: Proceedings of the SPIE 2563 Advanced Signal Processing Algorithms, San Diego, July 1995, pp. 2–20. [3] D.R. Brillinger, Time Series: Data Analysis and Theory, Holden-Day, San Francisco, CA, 1981. ¨ [4] P.-J. Chung, M.L. Jost, J.F. Bohme, Seismic wave parameter estimation and signal detection using broadband maximum likelihood methods, Computers and Geosciences 27 (2) (March 2001) 147–156. [5] P.-J. Chung, Robust ML estimation for unknown numbers of signals, in: Proceedings of the EUSIPCO, Poznan´, Poland, September 2007. ¨ ¨ [6] P.-J. Chung, J.F. Bohme, C.F. Mecklenbrauker, A.O. Hero, Detection of the number of signals using the Benjamini–Hochberg procedure, IEEE Transactions on Signal Processing 55 (6) (June 2007) 2497–2508. [7] P.-J. Chung, Stochastic maximum likelihood estimation under misspecified numbers of signals, IEEE Transactions on Signal Processing 55 (9) (September 2007) 4726–4731. ¨ [8] P.-J. Chung, C.F. Mecklenbrauker, Deterministic ML estimation for unknown numbers of signals, in: Proceedings of the EUSIPCO 2008, Lausanne, Switzerland, August 25–29, 2008. ¨ [9] P.-J. Chung, M. Viberg, C.F. Mecklenbrauker, Broadband ML estimation under model order uncertainty, in: Proceedings of the IEEE ICASSP 2009, Taipei, Taiwan, April 19–24, 2009. [10] P. Hall, The Bootstrap and Edgeworth Expansion, Springer, New York, 1992. [11] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985. [12] N. Johnson, S. Kotz, Continuous univariate distributions—2, in: Applied Probability and Statistics, Wiley, New York, 1970. ¨ [13] D. Kraus, Approximative Maximum–Likelihood–Schatzungen und ¨ verwandte Verfahren zur Ortung und Signalschatzung mit Sensorgruppen, Dr.-Ing. Dissertation, Ruhr University Bochum, Germany, Shaker Verlag, 1993. [14] D. Maiwald, Breitbandverfahren zur Signalentdeckung und -ortung mit Sensorgruppen in Seismik- und Sonaranwendungen, Dr.-Ing. Dissertation, Department of Electrical Engineering, Ruhr-Uni¨ Bochum, Germany, Shaker Verlag, 1995. versitat ¨ ¨ [15] C.F. Mecklenbrauker, P. Gerstoft, J.F. Bohme, P.-J. Chung, Hypothesis testing for geoacoustic environmental model using likelihood ratio, Journal of the Acoustical Society of America 105 (3) (March 1999) 1738–1748. [16] M. Wax, I. Ziskind, Detection of the number of coherent signals by the MDL Principle, IEEE Transactions on Acoustics, Speech, & Signal Processing 37 (8) (August 1989) 1190–1196. [17] M. Wax, I. Ziskind, Detection of the number of coherent signals by the MDL Principle, IEEE Transactions on Acoustics, Speech, & Signal Processing 37 (8) (August 1989) 1190–1196. [18] H. White, Maximum likelihood estimation under misspecified models, Econometrica 50 (1) (January 1982) 1–25. [19] K.M. Wong, Q.T. Zou, J.P. Reilly, P.C. Yip, On information theoretic criteria for determining the number of signals in high resolution array processing, IEEE Transactions on Acoustics, Speech, & Signal Processing 38 (11) (November 1990) 1959–1971.
Appendix. Cumulants of T ¼ logð1þðn1 =n2 ÞVÞ In this section, we derive the cumulants of the statistic T ¼ logð1 þ ðn1 =n2 ÞVÞ where V is Fn1 ;n2 distributed. Its inverse W ¼ 1=ð1þ ðn1 =n2 ÞVÞ ¼ w2n2 =w2n1 þ n2 is beta-distributed with parameters p ¼ n1 =2; q ¼ n2 =2 [12]. The pdf of W is given by fW ðwÞ ¼
1 wp1 ð1wÞq1 Bðp; qÞ
ð0 r wr 1Þ;
ð26Þ
where Bðp; qÞ represents the beta function Z 1 GðpÞGðqÞ Bðp; qÞ :¼ ¼ Bðq; pÞ: ð27Þ wp1 ð1wÞq1 dw ¼ Gðp þ qÞ 0 pffiffiffiffiffiffiffi Define j ¼ 1. The following expression for the characteristic function fT ðsÞ is valid if Refpg 4Imfsg, Refpg 40 and Refqg 4 0, Z 1 1 fT ðsÞ ¼ E½ejsT ¼ E½W js ¼ wp1 ð1wÞq1 wjs dw 0 Bðp; qÞ Z 1 1 Bðpjs; qÞ wq1js ð1wÞp1 dw ¼ : ð28Þ ¼ Bðp; qÞ 0 Bðp; qÞ Since p; q Z 12, s ¼ 0 belongs to the inner convergence region. By definition, the r th cumulant kr of T can be computed as follows: dr kr ¼ ðjÞr r logfT ðsÞ : ð29Þ ds s¼0 Combining (27) and (28), the logarithm of the characteristic function fT ðsÞ can be expressed as follows: logfT ðsÞ ¼ logBðq; pÞ þ logBðqjs; pÞ GðqjsÞGðpÞ ¼ logBðq; pÞ þlog Gðqjsþ pÞ ¼ logBðq; pÞ þ logGðpÞ þ logGðqjsÞlogGðqjs þ pÞ: ð30Þ
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Broadband ML estimation under model order uncertainty c ¨ Pei-Jung Chung a,, Mats Viberg b, Christoph F. Mecklenbrauker a
Institute for Digital Communications, The University of Edinburgh, UK Department of Signals and Systems, Chalmers University of Technology, Sweden c Institute of Communications and RF Engineering, Vienna University of Technology, Austria b
a r t i c l e i n f o
abstract
Article history: Received 8 June 2009 Received in revised form 5 November 2009 Accepted 6 November 2009 ¨ Dedicated to J.F. Bohme on the occasion of his 70th birthday Available online 26 November 2009
The number of signals hidden in data plays a crucial role in array processing. When this information is not available, conventional approaches apply information theoretic criteria or multiple hypothesis tests to simultaneously estimate model order and parameter. These methods are usually computationally intensive, since ML estimates are required for a hierarchy of nested models. In this contribution, we propose a computationally efficient solution to avoid this full search procedure and address issues unique to the broadband case. Our max-search approach computes ML estimates only for the maximally hypothesized number of signals, and selects relevant components through hypothesis testing. Furthermore, we introduce a criterion based on the rank of the steering matrix to reduce indistinguishable components caused by overparameterization. Numerical experiments show that despite model order uncertainty, the proposed method achieves comparable estimation and detection accuracy as standard methods, but at much lower computational expense. & 2009 Elsevier B.V. All rights reserved.
Keywords: Direction of arrival Broadband signals Maximum likelihood estimation Overparameterized models Unknown number of signals
1. Introduction The problem of estimating direction of arrival (DOA) plays an important role in array processing and a variety of applications such as radar, sonar, wireless communications and geophysics. The maximum likelihood (ML) approach is characterized by excellent statistical properties and robustness against small sample numbers, signal coherence and closely located sources. The standard ML approach assumes the number of signals, m, to be known and maximizes the likelihood function over an m-dimensional parameter set Mm . When m is unknown, traditional methods jointly estimate the number of signals and parameters of interest. The popular information theoretic criterion based approach [16,19] and the multiple hypothesis testing procedure [2,4,6,15]
Corresponding author.
E-mail addresses: [email protected] (P.-J. Chung), [email protected] (M. Viberg), [email protected] ¨ (C.F. Mecklenbrauker). 0165-1684/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.11.013
all belong to this category. More specifically, these methods compute ML estimates for a series of nested models and select the one that best fits the underlying criteria. Since ML estimates for all candidate models are required, the computational complexity becomes very high if the number of candidate models is large. Therefore, avoiding this full-search procedure can be a crucial step to improve computational efficiency. In this work, we pursue a max-search direction finding algorithm for unknown numbers of signals. The main idea is to consider the maximally hypothesized model and select relevant components associated with true parameters through simple hypothesis tests. In the narrow band case, the test can be conducted with an F-distributed statistic [8]. For broadband signals, the distribution under the null hypothesis has no closed form expression. Therefore, we apply the Cornish–Fisher expansion [10] to approximate the test threshold. In the previous publication [9], we exploited the important property that the overparameterized ML estimate contains the true parameters. This property was verified for narrow band, stochastic signal models [6]. The proof for the broadband,
ARTICLE IN PRESS P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
deterministic case will be presented here. Moreover, we shall discuss the identification problem caused by overparameterization and introduce a criterion to reduce indistinguishable components. The proposed procedure is different from the aforementioned model order estimation methods [6,14,17,19] in that it is aimed at extracting useful information about DOA parameters regardless of whether the number of signals is correctly specified or not. However, as a byproduct, the number of relevant components can be considered as an estimate for the actual number of signals. As it will be shown later, the proposed maxsearch procedure is computationally more attractive than the full search methods. The robustness against model order uncertainty could be of great importance to channel identification in wireless communications. In a time-varying propagation environment, our algorithm may help to distinguish the relevant point-like multipath components from the spurious ones with a minimum prior knowledge on model order. In the following, we give a brief description of the broadband signal model. Section 3 discusses asymptotic properties of ML estimation with overparameterized models. In Section 4, we develop the broadband ML estimation procedure and derive the Cornish–Fisher approximation to the test threshold. Simulation results are presented in Section 5, whereas Section 6 concludes the paper.
2. Broadband signal model Consider an array of n sensors receiving m broadband signals emitted by far-field sources located at h ¼ ½y1 ; . . . ; ym T . The sensor array outputs xðkÞ ðtÞ; ðt ¼ 0; . . . ; T1Þ are divided into K snapshots of length T 0 ¼ T=K and short-time Fourier-transformed T 1 1 X X ðkÞ ðoÞ ¼ pffiffiffi wðtÞxðkÞ ðtÞejot ; Tt¼0
ð1Þ
where fwðtÞgT1 t ¼ 0 is a window function. For large sample size, T 0 -1, the frequency domain data are approximately described by the following relation: X ðkÞ ðoÞ ¼ H m ðo; hm ÞS ðkÞ ðoÞ þ U ðkÞ ðoÞ;
ð2Þ
where the matrix H m ðo; hm Þ ¼ ½dðo; y1 Þ dðo; yi Þ dðo; ym Þ 2 Cnm consists of m steering vectors with the i th column dðo; yi Þ, corresponding to the i th incoming wave. We assume that the array manifold is unambiguous, i.e. any collection of m steering vectors corresponding to distinct parameters is linearly independent, and smooth, i.e. the third derivatives of dðo; yi Þ with respect to yi are bounded. The signal waveform S ðkÞ ðoÞ is considered as unknown and deterministic. The spatially white noise U ðkÞ ðoÞ results only from sensors. Under certain regularity conditions [3], X ðkÞ ðoj Þ ðk ¼ 1; . . . ; K; j ¼ 1; . . . ; JÞ are asymptotically independent, complex normally distributed with mean H m ðoj ; hm ÞS ðkÞ ðoj Þ and covariance matrix nðoj ÞI,
1351
where nðoj Þ is the unknown noise spectral parameter and I is an identity matrix of conformable dimension. Given a pre-specified number of signals, m, the broadband ML estimate is obtained by minimizing the negative log-likelihood function over the m-dimensional parameter space Mm [1,13]:
h^ m ¼ arg min lðhm Þ; hm 2Mm
lðhm Þ ¼
J X
logtr½ðIPðoj ; hm ÞÞC^ x ðoj Þ;
ð3Þ
j¼1
where Pðoj ; hm Þ represents the projection matrix onto the subspace spanned by the columns of H m ðoj ; hm Þ. The power spectral estimate of array outputs is give by P C^ x ðoj Þ ¼ ð1=KÞ Kk ¼ 1 X ðkÞ ðoj ÞX ðkÞ ðoj ÞH . In many array processing applications, the number of signals m is unknown. The most successful approaches exploit the information theoretic criteria or carry out multiple hypothesis testing to jointly estimate the DOA parameters and the model order. Given a hierarchy of nested models represented by nested parameter spaces Mm , i.e. hm 2 Mm for 1 r m rM, M 1 M2 MM ;
ð4Þ
these methods derive ML estimates for all candidate models and select the best one according to the underlying decision rule. Such full search procedures are typically computationally demanding, due to the multidimensional nonlinear minimization required by each candidate model. In contrary to this full search idea, we suggest a computationally more efficient alternative that considers the maximally hypothesized model MM and extracts DOA information from the corresponding estimate. 3. ML estimation of overparameterized models The key idea behind the proposed algorithm is that asymptotically, for an overparameterized model, the ML estimate contains components of the true DOA parameter, ho . By overparameterization, we mean that the assumed number of signals m is larger than the true number of signals mo . In other words, the asymptotic estimate derived from m 4mo ; hm ¼ ½y1 ; y2 ; . . . ; ym T , contains mo components associated with those of ho ¼ ½yo1 ; yo2 ; . . . ; yomo T . The remaining ðmmo Þ components correspond to redundant parameters. Applying the asymptotic theory of quasi-ML estimation [18], the analysis in [7] shows that under stochastic signal models, the narrow band ML estimator converges to the minimizing point of the Kullback–Leibler distance between the pdf of the assumed model and that of the true model. This overparameterized estimate hm is characterized by the following asymptotic property (valid for large K): spðH m ðo; hm ÞÞ*spðH mo ðo; ho ÞÞ;
ð5Þ
where spðH m ðo; hm ÞÞ and spðH mo ðo; ho ÞÞ denote the signal subspaces corresponding to hm and ho , respectively.
ARTICLE IN PRESS 1352
P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
Clearly, hm with the true parameter ho as a subvector is a solution to (5). These results can be generalized to broadband ML estimation with the deterministic (or conditional) signal model. However, rather than minimizing the Kullback– Leibler distance as in the stochastic case, we shall minimize the asymptotic likelihood function directly. This approach does not require standard assumptions on ML estimation such as fixed number of unknown parameters and i.i.d. samples which are not fulfilled in the current signal model. To investigate large sample properties, we shall exploit the true covariance matrix of array outputs H
C xo ðoÞ ¼ H mo ðo; ho ÞC so ðoÞH mo ðo; ho Þ þ no ðoÞI;
ð6Þ
where C so ðoÞ ¼ E½S ðkÞ ðoÞS ðkÞ ðoÞH and no ðoÞ represent the signal covariance matrix and noise level, respectively. It is well established that for K-1, the sample covariance matrix C^ x ðoÞ converges with probability one to the ensemble covariance matrix C xo ðoÞ. Result 1. Let m denote the assumed number of signals and m4 mo . For large sample size, the ML estimate obtained by minimizing the broadband likelihood function (3) converges to hm that minimizes the criterion Q ðhm Þ ¼
J X
logtr½ðIPðoj ; hm ÞÞC xo ðoj Þ
ð7Þ
j¼1
and achieves the minimum Q ðhm Þ ¼
J X
lognðoj Þ þ JlogðnmÞ:
ð8Þ
j¼1 j ; hm ÞÞ
Furthermore, the assumed signal subspace spðH m ðo and the true signal subspace spðH mo ðoj ; ho ÞÞ are asymptotically related as follows (valid for large K): spðH m ðoj ; hm ÞÞ*spðH mo ðoj ; ho ÞÞ;
j ¼ 1; . . . ; J:
ð9Þ
Proof. Recall that the finite sample likelihood function (3) is given by lðhm Þ ¼
J X
logtr½ðIPðoj ; hm ÞÞC^ x ðoj Þ:
rewritten as Q ðhm Þ ¼
J X
Q ðhm ; oj Þ;
j¼1
Q ðhm ; oj Þ ¼ logtr ½U 2 ðoj ÞH C xo ðoj ÞU 2 ðoj Þ:
ð10Þ
According to Corollary 4.3.18 of [11], the minimum of Q ðhm ; oj Þ subject to the constraint U 2 ðoj ÞH U 2 ðoj Þ ¼ I nm Q ðhm ; oj Þ ¼ min logtr ½U 2 ðoj ÞH C xo ðoj ÞU 2 ðoj Þ U 2 ðoj Þ
¼ log½ðnmÞnðoj Þ
ð11Þ
is achieved if the columns of the minimizing matrix U 2 ðoj Þ are orthonormal eigenvectors corresponding to the ðnmÞ smallest eigenvalues of C xo ðoj Þ. The ðnmÞ smallest eigenvalues are equal to the noise spectral parameter nðoj Þ because ðnmÞ oðnmo Þ. This also implies that spðU 2 ðoj ÞÞ is a subspace of the noise subspace of C xo ðoj Þ. Since U ðoj Þ ¼ ½U 1 ðoj Þ U 2 ðoj Þ forms an orthonormal basis for Cn1 , spðU 1 ðoj ÞÞ*spðH mo ðoj ; ho ÞÞ:
ð12Þ
The relation (9) is a direct consequence of spðU 1 ðoj ÞÞ ¼ spðHðoj ; hm ÞÞ and it is satisfied for j ¼ 1; . . . ; J if hm includes the true parameter ho as a subvector. As all frequency bins are minimized at hm , the broadband criterion Q ðhm Þ achieves its minimum (8) at hm . The minimum asymptotic likelihood is given by Q ðhm Þ ¼
J X
lognðoj Þ þ JlogðnmÞ:
j¼1
This completes the proof.
&
Eq. (8) shows the minimum asymptotic likelihood Q ðhm Þ varies with the assumed number of m. In contrast to the stochastic model whose asymptotic likelihood remains constant for m 4mo [7], the selection of relevant estimates can no longer rely on each component’s contribution to the likelihood function. Therefore, we propose a solution based on hypothesis testing in the following section.
j¼1
Following the fact C^ x ðoj Þ-C xo ðoj Þ for K-1, the asymptotic likelihood function becomes Eq. (7). Provided the array response (steering vector) is a smooth function of hm , it is straightforward to verify that the convergence is uniform in the parameters. Consequently, the parameter estimates converge to the minimizing arguments of (7) in the same statistical sense (e.g. w.p. 1) as the sample covariance matrix converges. To find the minimum of Q ðhm Þ, consider the singular value decomposition H m ðoj Þ ¼ Uðoj ÞRðoj ÞVðoj ÞH where Uðoj Þ ¼ ½U 1 ðoj Þ U 2 ðoj Þ is a unitary matrix. The n m matrix U 1 ðoj Þ consists of the first m columns of Uðoj Þ corresponding to the m largest singular values. Then we can express the projection matrix as Pðoj Þ ¼ U 1 ðoj Þ U 1 ðoj ÞH and its orthogonal complement as IPðoj Þ ¼ U 2 ðoj ÞU 2 ðoj ÞH . The asymptotic likelihood (7) can be
4. The max-search approach Motivated by the above observation, we suggest to compute the ML estimate for the maximally hypothesized model, m ¼ M, and select relevant components as estimates. More precisely, the proposed algorithm considers the estimate derived from an M-dimensional space
h^ M ¼ argminlðhM Þ; hM
ð13Þ
where M Zmo . The M 1 vector h^ M ¼ ½y^ 1 ; . . . ; y^ M T contains more elements than the mo 1 vector ho . The elements of h^ M that are associated with those of ho , are referred to as relevant components. The remaining ðMmo Þ components of h^ M are the redundant components. Clearly, the key step in the proposed max-search procedure is the identification of relevant components.
ARTICLE IN PRESS P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
In the stochastic model, relevant components are selected by thresholding the likelihood function [5], because redundant components do not contribute to the likelihood function. For the deterministic case, the asymptotic minimum (8) varies with assumed model order m and does not allow such simple selection scheme. We suggest the following hypothesis test to validate the ith component:
AM : XðoÞ ¼ HM ðo; h^ M ÞS M ðoÞ þUðoÞ;
ð15Þ
j¼1
J 1X n1 log 1 þ F i ð oj Þ : n2 J j¼1
ð16Þ
Here, Fi ðoj Þ can be interpreted as an incremental signalto-noise ratio F i ð oj Þ ¼
n2 tr½ðPðoj ; h^ M ÞPðoj ; h~ i ÞÞC^ x ðoj Þ n1 tr½ðIPðoj ; h^ M ÞÞC^ x ðoj Þ
and it turns out to be asymptotically Fn1 ;n2 distributed (valid for large K) under suitable assumptions (see Section 4.1). The associated degrees of freedom are [14] n1 ¼ 3K;
n2 ¼ Kð2n2M1Þ:
ð17Þ
The component y^ i is relevant if Hi is rejected. Given a significance level a, the null hypothesis Hi is rejected if Ti exceeds a threshold ta . Consequently,
y^ i is relevant if Ti Z ta :
ð18Þ
The output of the algorithm is the relevant vector that contains all components
h^ o ¼ ½y^ ð1Þ ; . . . ; y^ ðpÞ T :
(2) Compute test statistic Ti , i ¼ 1; . . . ; M. (3) Test Hi against AM , i ¼ 1; . . . ; M. y^ is relevant if T Z ta . i
ð14Þ
contains all elements of h^ M except the ith component. The n ðM1Þ matrix H M1 ðh~ i Þ contains steering vectors corresponding to the DOA parameters in h^ M1 . Given the estimate h^ M , the hypothesis test (14) decides whether the signal Si ðoÞ associated with the ith component is zero. Note that Si ðoÞ ¼ 0 implies y^ i does not correspond to any actual signal source. The test statistic Ti ði ¼ 1; 2; . . . ; MÞ derived from the likelihood ratio principle, lT ðh~ i ÞlT ðh^ M Þ, is given by ! J tr½ðIPðoj ; h~ i ÞÞC^ x ðoj Þ 1X Ti ¼ log J tr½ðIPðoj ; h^ M ÞÞC^ x ðoj Þ ¼
Input: fX ðkÞ ðoj Þ: k ¼ 1; . . . ; K; j ¼ 1; . . . ; Jg, upper bound on number of signals: M, test threshold: ta . (1) Find the ML estimate h^ M .
Output: relevant vector h^ o ¼ ½y^ ð1Þ ; y^ ð2Þ ; . . . ; y^ ðpÞ T .
where Hi and AM represent the null hypothesis and the alternative, respectively. The ðM1Þ 1 vector
h~ i ¼ ½y^ 1 y^ i1 y^ i þ 1 y^ M T
Table 1 The max-search algorithm for broadband ML estimation with unknown numbers of signals.
i
Hi : XðoÞ ¼ HM1 ðo; h~ i ÞS~ M1 ðoÞ þ UðoÞ;
1353
ð19Þ
As a byproduct of the algorithm, the number of relevant components provides an estimate for the number of signals. However, we emphasize that the primary concern of the proposed approach is parameter estimation. Unlike methods designed for model order determination, whose performance is measured solely by correctness of the number of signals, the output of the proposed algorithm h^ o still contains useful information about the true parameters ho even when the number of signals is misspecified. Given an upper bound on the
number of the signals M and the threshold ta , the proposed max-search algorithm is summarized in Table 1.
4.1. Cornish–Fisher approximation Under the null hypothesis Hi , the statistic Fi ðoj Þ is asymptotically Fn1 ;n2 distributed with the degrees of freedom n1 , n2 given by (17) for large K. As (20) shows, in the narrow band case, J ¼ 1, the test can be equivalently conducted with an F-statistic. However, for broadband signals, the distribution under the null hypothesis has no analytical expression. To overcome this difficulty, we suggest the Cornish–Fisher expansion to approximate the test threshold. The summands in (16) can be interpreted as i.i.d. samples drawn from the random variable n1 ð20Þ Fi ; T i ¼ log 1 þ n2 where the statistic F i is Fn1 ;n2 - distributed. This property follows from the asymptotic independence of the Fouriertransformed data X ðkÞ ðoj Þ ðj ¼ 1; . . . ; JÞ and the same degrees of freedom, n1 and n2 , for each frequency bin. Let T i ¼ J 1=2
ðT i mT Þ
sT
ð21Þ
denote the standardized statistic where mT ¼ E½T i and s2T ¼ E½T i mT 2 are the mean and variance of T i , respectively. With the knowledge of mean and variance, the distribution of T i can be approximated by the Edgeworth expansion [10]. For a statistic whose distribution admits an Edgeworth expansion, the a- quantile, defined by PfT i rt a g ¼ 1a
ð22Þ
can be approximated by the Cornish–Fisher expansion, 1 1 1 t a ¼ za þ pffiffi p1 ðza Þ þ p2 ðza Þ þ þ r=2 pr ðza Þ þ oðJr=2 Þ; J J J ð23Þ where za is the standard normal aquantile. The polynomial pr ðzÞ is of degree at most r þ 1, and depend on the cumulants only up to order r þ2. This expansion is valid for fixed r, as J-1. The quality of the approximation is given by the remainder oðJr=2 Þ. For fixed samples, the accuracy of the approximation can be improved by including more terms in the expansion. The first two
ARTICLE IN PRESS 1354
P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
polynomials are given by p1 ðzÞ ¼
k3 6
Sample mean of relevant components 37
ðz2 1Þ;
36
k 3 2 k 4 k4 2 z ðz 1Þ þ z ðz 3Þ þ ðz 10z þ 15Þ : p2 ðzÞ ¼ 18 24 72 2 3
35
r
where Cðr1Þ ðsÞ ¼ dr =ds ðlogGðsÞÞ denotes the r th derivative of the logarithm of the gamma function. The mean mT and variance s2T are related to the cumulants as follows:
s2T ¼ k2 :
33 32 31 30 29
ð25Þ
28
Computational details of cumulants kr are presented in Appendix.
27
mT ¼ k1 ;
θ(1), M = 4 θ(2), M = 4 θ(1), M = 3 θ(2), M = 3
34
The r th cumulant kr of T i depends on the degrees of freedom n1 , n2 through the following formula: h n n n i kr ¼ ð1Þr Cðr1Þ 1 Cðr1Þ 1 þ 2 ; ð24Þ 2 2 2
Mean (degree)
2 3
−10 −8
−6
−4
−2
0
2
4
6
8
10
SNR (dB)
4.2. Identification of indistinguishable components One situation that may occur in the proposed algorithm is that the estimate h^ M contains indistinguishable components, which leads to a rank deficient steering matrix Hðo; h^ M Þ. By ‘‘indistinguishable components’’ we mean that two DOA parameters are too close or almost identical at the minimum of the likelihood function. Under the assumption of unambiguous array manifold, property (9) permits such parameters because the underlying model is overparameterized. In this case, both the relevant vector h^ o and the test (14) no longer provide satisfying results. To overcome difficulties caused by overparameterization, we suggest to examine the smallest singular value s^ M of the steering matrix Hðo; h^ M Þ.1 If js^ M j r d where d is a pre-specified small positive number, then the procedure computes the ML estimate for the next reduced model MM1 . If the resulting Hðo; h^ M1 Þ is full rank, then the relevant components are chosen from h^ M1 . Otherwise, the same step is repeated for the next reduced model. In the simulation, we observe that this modified step is carried out only by 5% of all trials. The computational cost is still much lower than the full search approach. The indistinguishable components are usually almost identical; for example, the second and third components of the h^ 4 ¼ ½37:93783 29:21293 29:21353 60:39263 estimate differ from each other by less than 0:0063 . Hence, it is sufficient to consider the steering matrix at one frequency bin. In Section 5, the threshold d is chosen based on numerical data. How to select d in a statistically justified manner is still under investigation. 5. Simulation In the simulation, a uniform linear array of 10 sensors with inter-element spacings of half a wavelength l=2 is employed. The wavelength is given by l ¼ v=f0 where v represents the propagation velocity and f0 is a pre-selected 1 Alternatively, a minimum separation in the elements of h^ M could be used.
Fig. 1. Sample mean of relevant components. M ¼ 3, 4. The true DOA parameter ho ¼ ½283 363 , SNR ¼ ½10 : 2 : 10 dB.
reference frequency. The broadband signals are generated by m0 ¼ 2 uncorrelated signals, located at h0 ¼ ½283 363 with respect to array broadside, and of various strengths. We choose J ¼ 11 frequency bins equally spaced between ð17=32Þf0 and f0 for processing. The number of snapshots is K ¼ 50. The difference in signal strengths is [1 0] dB, where 0 dB corresponds to the reference signal. The signal to noise ratio (SNR), defined as 10 log10 E½jSi ðoj Þj2 =nðoj Þ for the i th signal ði ¼ 1; 2Þ, varies from 10 to 10 dB in a 2 dB step. We consider two upper bounds on the number of signals: M ¼ 3 and 4. The latter represents a larger mismatch in the model order. Each experiment consists of 200 Monte Carlo trials. The test level is chosen to be a ¼ 0:05. The parameter d ¼ 0:02 is used to test whether the steering matrix is rank deficient or not. Fig. 1 shows the sample mean of the relevant components y^ ð1Þ , y^ ð2Þ , respectively. For both M ¼ 3 and 4, the bias is less than 0.061 over the entire SNR range. Since the results are obtained from finite samples, we conjecture that y^ ð1Þ and y^ ð2Þ are asymptotically bias free. The empirical variance of y^ ð1Þ and y^ ð2Þ is presented in Fig. 2. For comparison, we also show results obtained from the no-mismatch case with M ¼ m0 ¼ 2. All three curves decline with increasing SNR. The estimates obtained from M ¼ m0 ¼ 2 have the smallest variance. For both relevant components y^ ð1Þ and y^ ð2Þ , M ¼ 4 results in a larger variance than M ¼ 3. This suggests that the variance increases with a larger degree of mismatch. In addition, it is observed in the simulation that bias and variance can be reduced by increasing the number of data samples K or frequency bins J. In Fig. 3, we compare the probability of correct detection of the proposed algorithm with the multiple hypothesis testing procedure [6]. By ‘‘correct detection’’, we mean that the number of relevant components equals the true number of signals. At the low SNR region, 10 to 4 dB, the curves below 100% indicate that the signals are not always detected. However, both M ¼ 3 and 4 have a higher probability of correct detection than the multiple test. For SNR 6 to 10 dB,
ARTICLE IN PRESS P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
Empirical variance of θ(1)
101 Variance (degree2)
1355
100
θ(1), M = 4
10−1
θ(1), M = 3 θ(1), m0 = 2
10−2 −10
−8
−6
−4
0 SNR (dB)
2
4
6
8
10
4
6
8
10
Empirical variance of θ(2)
101 Variance (degree2)
−2
100 θ(2), M = 4
10−1
θ(2), M = 3 θ(2), m0 = 2
10−2 −10
−8
−6
−4
−2
0
2
SNR (dB) Fig. 2. Empirical variance of relevant estimates y^ ðiÞ , i ¼ 1, 2.
signals is unknown or misspecified. The price paid for model uncertainty is increased variance. Furthermore, the comparison between the proposed algorithm and the full search procedure [6] showed that without computing ML estimates for each candidate model, the former has a lower SNR threshold and achieves good detection performance at high SNRs. The computational time required by the suggested max-search procedure is on average 35% less than the full search procedure. The saving is expected to grow for larger numbers of candidate models.
Estimation of Number of Signals 1 Probability of Correct Detection
0.9 0.8
M=4 M=3 Multiple test
0.7 0.6 0.5 0.4 0.3
6. Conclusion
0.2 0.1 0 −10 −8
−6
−4
−2
0
2
4
6
8
10
SNR (dB) Fig. 3. Probability of correct detection. Comparison of multiple test procedure [5], M ¼ 3, 4.
the multiple test procedure achieves 100% probability of correct detection, while M ¼ 3; 4 increase from 90% to 95%. The curve associated with M ¼ 3 shows a slightly higher probability of correct detection than that of M ¼ 4. In summary, simulation results demonstrate that relevant components convey useful information about the true parameters even when the correct number of
We developed a broadband ML estimation procedure for unknown numbers of signals. The suggested algorithm computes ML estimates only for the maximally hypothesized number of signals. The relevant components associated with the true DOA parameters are selected by simple hypothesis tests. The test threshold is approximated by the Cornish–Fisher expansion. We also introduce a criterion to reduce indistinguishable components caused by overparameterization. Compared to traditional methods for joint parameter estimation and signal detection, the proposed max-search approach avoids the full search process through a series of nested models, which leads to significant improvement in computational efficiency if the number of candidate models is large. Numerical results showed that the proposed algorithm achieves comparable estimation accuracy
ARTICLE IN PRESS 1356
P.-J. Chung et al. / Signal Processing 90 (2010) 1350–1356
as the standard ML approach does. The number of signals can also be accurately determined by the number of relevant components. The proposed algorithm provides a computationally attractive alternative to existing methods, particularly in the low SNR region. A potential drawback of the proposed approach is that over-specification of the number of signals may lead to higher sensitivity to errors in the steering vector model. This is an interesting topic for further investigations.
The r th derivative of (30) leads to the desired r th cumulant kr ðr ¼ 1; 2; . . .Þ of T dr kr ¼ r logfT ðsÞ ds s¼0 r d dr ¼ log G ðqjsÞ log G ðqjs þ pÞ r r ds ds s¼0 ¼ ð1Þr ½Cðr1Þ ðqÞCðr1Þ ðq þ pÞ:
ð31Þ
We obtain (24) by replacing p ¼ n1 =2, q ¼ n2 =2 into the above equation.
Acknowledgements
References
The authors wish to thank the anonymous reviewers for their comments that significantly improved the manuscript and also thank the Managing Editor Prof. A.M. Zoubir for coordinating a speedy review.
¨ [1] J.F. Bohme, Array processing, in: S. Haykin (Ed.), Advances in Spectrum Analysis and Array Processing, Prentice-Hall, Englewood Cliffs, NJ, 1991, pp. 1–63. ¨ [2] J. F. Bohme, Statistical array processing of measured sonar and seismic data, in: Proceedings of the SPIE 2563 Advanced Signal Processing Algorithms, San Diego, July 1995, pp. 2–20. [3] D.R. Brillinger, Time Series: Data Analysis and Theory, Holden-Day, San Francisco, CA, 1981. ¨ [4] P.-J. Chung, M.L. Jost, J.F. Bohme, Seismic wave parameter estimation and signal detection using broadband maximum likelihood methods, Computers and Geosciences 27 (2) (March 2001) 147–156. [5] P.-J. Chung, Robust ML estimation for unknown numbers of signals, in: Proceedings of the EUSIPCO, Poznan´, Poland, September 2007. ¨ ¨ [6] P.-J. Chung, J.F. Bohme, C.F. Mecklenbrauker, A.O. Hero, Detection of the number of signals using the Benjamini–Hochberg procedure, IEEE Transactions on Signal Processing 55 (6) (June 2007) 2497–2508. [7] P.-J. Chung, Stochastic maximum likelihood estimation under misspecified numbers of signals, IEEE Transactions on Signal Processing 55 (9) (September 2007) 4726–4731. ¨ [8] P.-J. Chung, C.F. Mecklenbrauker, Deterministic ML estimation for unknown numbers of signals, in: Proceedings of the EUSIPCO 2008, Lausanne, Switzerland, August 25–29, 2008. ¨ [9] P.-J. Chung, M. Viberg, C.F. Mecklenbrauker, Broadband ML estimation under model order uncertainty, in: Proceedings of the IEEE ICASSP 2009, Taipei, Taiwan, April 19–24, 2009. [10] P. Hall, The Bootstrap and Edgeworth Expansion, Springer, New York, 1992. [11] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985. [12] N. Johnson, S. Kotz, Continuous univariate distributions—2, in: Applied Probability and Statistics, Wiley, New York, 1970. ¨ [13] D. Kraus, Approximative Maximum–Likelihood–Schatzungen und ¨ verwandte Verfahren zur Ortung und Signalschatzung mit Sensorgruppen, Dr.-Ing. Dissertation, Ruhr University Bochum, Germany, Shaker Verlag, 1993. [14] D. Maiwald, Breitbandverfahren zur Signalentdeckung und -ortung mit Sensorgruppen in Seismik- und Sonaranwendungen, Dr.-Ing. Dissertation, Department of Electrical Engineering, Ruhr-Uni¨ Bochum, Germany, Shaker Verlag, 1995. versitat ¨ ¨ [15] C.F. Mecklenbrauker, P. Gerstoft, J.F. Bohme, P.-J. Chung, Hypothesis testing for geoacoustic environmental model using likelihood ratio, Journal of the Acoustical Society of America 105 (3) (March 1999) 1738–1748. [16] M. Wax, I. Ziskind, Detection of the number of coherent signals by the MDL Principle, IEEE Transactions on Acoustics, Speech, & Signal Processing 37 (8) (August 1989) 1190–1196. [17] M. Wax, I. Ziskind, Detection of the number of coherent signals by the MDL Principle, IEEE Transactions on Acoustics, Speech, & Signal Processing 37 (8) (August 1989) 1190–1196. [18] H. White, Maximum likelihood estimation under misspecified models, Econometrica 50 (1) (January 1982) 1–25. [19] K.M. Wong, Q.T. Zou, J.P. Reilly, P.C. Yip, On information theoretic criteria for determining the number of signals in high resolution array processing, IEEE Transactions on Acoustics, Speech, & Signal Processing 38 (11) (November 1990) 1959–1971.
Appendix. Cumulants of T ¼ logð1þðn1 =n2 ÞVÞ In this section, we derive the cumulants of the statistic T ¼ logð1 þ ðn1 =n2 ÞVÞ where V is Fn1 ;n2 distributed. Its inverse W ¼ 1=ð1þ ðn1 =n2 ÞVÞ ¼ w2n2 =w2n1 þ n2 is beta-distributed with parameters p ¼ n1 =2; q ¼ n2 =2 [12]. The pdf of W is given by fW ðwÞ ¼
1 wp1 ð1wÞq1 Bðp; qÞ
ð0 r wr 1Þ;
ð26Þ
where Bðp; qÞ represents the beta function Z 1 GðpÞGðqÞ Bðp; qÞ :¼ ¼ Bðq; pÞ: ð27Þ wp1 ð1wÞq1 dw ¼ Gðp þ qÞ 0 pffiffiffiffiffiffiffi Define j ¼ 1. The following expression for the characteristic function fT ðsÞ is valid if Refpg 4Imfsg, Refpg 40 and Refqg 4 0, Z 1 1 fT ðsÞ ¼ E½ejsT ¼ E½W js ¼ wp1 ð1wÞq1 wjs dw 0 Bðp; qÞ Z 1 1 Bðpjs; qÞ wq1js ð1wÞp1 dw ¼ : ð28Þ ¼ Bðp; qÞ 0 Bðp; qÞ Since p; q Z 12, s ¼ 0 belongs to the inner convergence region. By definition, the r th cumulant kr of T can be computed as follows: dr kr ¼ ðjÞr r logfT ðsÞ : ð29Þ ds s¼0 Combining (27) and (28), the logarithm of the characteristic function fT ðsÞ can be expressed as follows: logfT ðsÞ ¼ logBðq; pÞ þ logBðqjs; pÞ GðqjsÞGðpÞ ¼ logBðq; pÞ þlog Gðqjsþ pÞ ¼ logBðq; pÞ þ logGðpÞ þ logGðqjsÞlogGðqjs þ pÞ: ð30Þ