Efficient wideband signals’ direction of arrival estimation method with unknown number of signals

Abstract

This article focuses on the direction of arrival estimation of wideband signals, and an efficient method based on test of orthogonality of frequency subspaces is proposed. Using suitable coefficient m to weight and inverse the covariance matrix, the proposed method can obtain estimated noise subspace with unknown number of signals. Unlike conventional subspace estimation with heavy complex computations, the proposed method reduces the amount of computation by real matrix multiplication, and through only half of the total angular field-of-view. It decreases the computational load of spectrum peak search to half. Moreover, compared with the existing real-valued techniques suitable for only the centro-symmetrical array, the new method can be applied to arbitrary array. The simulation and experimental data results are provided to demonstrate the effectiveness and performance of the proposed method.

Introduction

Source localization using wideband signals has been widely studied in multiple areas with very diverse application such as radar, communication, and microphone array systems,^1,2 and a variety of techniques have been proposed to handle the wideband signal localization problem. The maximum likelihood (ML) method is extended to the signals’ direction of arrival (DOA) estimation in Ye and Degroat ³ and Souza and Fortes; ⁴ since the ML method involves to maximize the likelihood function with highly nonlinearity, its computational complexity is very heavy. The most representative wideband DOA estimation approaches are incoherent signal subspace method (ISSM) ⁵ and coherent signal subspace method (CSSM).^6–8 ISSM estimates sources’ DOAs separately at each frequency bin and then joins these results to get a final estimate. However, it cannot deal with coherent signals and has bad performance at low signal-to-noise ratio (SNR). CSSM is the most typical coherent wideband focus algorithm and it converts the wideband signal subspace into a predefined narrowband subspace by focusing matrices through suitable transformation matrices, and subsequently the narrowband subspace-based DOA estimation methods, for example, MUSIC can be used. CSSM is effective in low SNR cases, while it is sensitive to the initial values and poor initial values can result in bigger biased estimates. The focusing methods can reduce the estimation bias; however, their performance extremely relies on the accuracy of initial pre-estimation angles. To overcome the problem of pre-estimation angle, Yoon et al. ⁹ proposed test of orthogonality of projected subspaces (TOPS) algorithm which does not require focusing angles, but its performance depends on the selection of reference frequency, and when the reference frequency is unsuitable, the performance of TOPS algorithm would become worse. Zhang et al. ¹⁰ presented extended TOPS method which has better estimation performance than TOPS; however, how to choose reference frequency is still an open problem. Yu et al. ¹¹ presented test of orthogonality of frequency subspaces (TOFS) method which exploits the information of multiple frequencies simultaneously and does not require any preprocessing for initial values. TOFS can eliminate the pseudo peak of TOPS and has better performance than TOPS method; however, it needs exhaustive spectral search which limits its practical application. In addition, there exist intrinsic limitations of the subspace decomposition approaches, that is, a priori information of incident signal number needs to be known in advance. The information theoretic criteria such as Akaike information criterion (AIC) and minimum description length (MDL) ¹² and their variants^13,14 are the most significant approaches for detecting the number of sources. For Qi et al., ¹³ it requires that the number of sources is no bigger than half of the number of sensors, and for Huang et al., ¹⁴ it is not available in the case of coherent signals. Qian et al. ¹⁵ proposed a MUSIC-like algorithm which does not need to know the number of signal; however, it has lower estimation accuracy because of the reconstruction of Toeplitz matrix. Reddy et al. ¹⁶ proposed a low-complexity algorithm with unknown number of sources but only limit with the use of narrowband signals. Recently, the emerging field of sparse representation of sources has attracted tremendous attention, and many researchers have shown enormous interests for its perfect resolution.^17–19 Via covariance matrix sparse representation, DOA estimation can be performed where the prior information of the signal number is also not needed. Using the idea of sparse representation, wideband DOA estimation has been realized in Xu et al. ²⁰ and Zhao et al.; ²¹ however, these sparse recovery algorithms need complex optimization process and spectral search. Generally, the computational complexity of spectral search is $O(J{M}^2)$ , where M is the sensor number and J stands for the total sample points for spectral search; usually, J is larger than M, so the heavy complexity of spectral search limits their application in practice.

Note that all the wideband signals’ localization methods aforementioned share a common problem, that is, their computational cost is very high. The computation of wideband DOA estimation methods based on subspace decomposition consists of two parts: one is the calculation of the estimated noise subspace and the other is the search of spectrum peak. Generally, the array covariance matrix is a complex matrix, since the complex multiplication costs four times more than that of real multiplications, and a unitary transform can convert the complex matrix into a real matrix along with their eigenvectors and thereby decreasing the amount of calculation of the estimated subspace at least by a factor of 4 without losing accuracy. ²² Haardt and Nossek ²³ applied this technique to Estimation of Signal Parameters via Rotation Invariance Technique (ESPRIT) to reduce the computational burden, and Yilmazer et al. ²⁴ extended this method into matrix pencil (MP) method to successfully decrease the computational cost. Yan et al. ²⁵ proposed the semi-real-valued Capon method which could reduce 50% computational complexity of spectrum search and thus, significantly reduce the computation of peak search.

In this article, we propose an efficient DOA estimation method for wideband signals with unknown number of signals, which can overcome the aforementioned shortcomings of the existing algorithm. First, the proposed algorithm can obtain the real sample covariance matrix of each frequency bin through the addition of covariance matrix and its conjugate, which can reduce the computational load of matrix multiplication by about 75%. Then, we can obtain the estimated noise subspace by reversing the weighted covariance matrix through suitable coefficient m with no prior information of the source number. Finally, using symmetrical characteristics of new spectrum, just through only half of the total angular field-of-view, the proposed method can acquire wideband signals’ DOA estimation very fast.

This article is organized as follows. Section “Wideband signals’ model” provides wideband signals’ model that will be used through the article. Section “Test of orthogonality of subspaces” introduces test of the orthogonality of the subspace theory where TOPS and TOFS are described. A description of the modified TOFS (MTOFS) method is given in section “MTOFS,” followed by the numerical simulation results and discussion in section “Simulation results.” Finally, conclusions are given in section “Conclusion.”

Notation: We use boldface upper-case letters (e.g. A ) to denote matrices and boldface lower-case letters (e.g. a ) to denote vectors. $E(•)$ represents the statistical mean operation. The superscripts $(•{)}^{*},(•{)}^T,(•{)}^H$ , and $(•{)}^{-1}$ represent conjugate, transpose, Hermitian transpose, and inverse, respectively. ${\mathit{I}}_M$ denotes the $M\times M$ identity matrix and the symbol $\mathrm{diag}\{z_1,z_2,\dots ,z_m\}$ presents a diagonal matrix with diagonal entries $z_1,z_2,\dots ,z_m$ .

Wideband signals’ model

Consider a scenario where P wideband coherent signals impinge on a uniform linear array (ULA) which is composed of M sensors. The bandwidths of the wideband signals need not be identical, but there should be some frequency bandwidth $[f_L,f_H]$ and the signals come from directions ${\theta }_p,p=1,\dots ,P$ . The $m\mathrm{th}$ sensor output is $x_m(t)$ which can be written as:

$$x_m(t)=\sum _{p=1}^P{s}_p[t-{\tau }_m({\theta }_p)]+n_m(t),1\le m\le M$$

(1)

where $s_p(t)$ is the pth wideband signal, $n_m(t)$ is the noise loading at the mth sensor which is assumed to be uncorrelated with the signals and white both temporally and spatially, and ${\tau }_m({\theta }_p)$ is the propagation delay associated with the mth sensor and pth source which is defined by

$${\tau }_m({\theta }_p)=\frac{(m-1)d\sin ({\theta }_p)}{c}$$

(2)

where c is the velocity of source propagation.

Suppose all sources could be partitioned into J non-overlapping narrowband blocks and the observation time is long enough so that the Fourier transform of the sensor output has good resolution. Via discrete-time Fourier transform (DTFT), the array output at frequency bin j can be expressed as

$$\mathit{X}(f_j)=\mathit{A}(f_j,\theta )\mathit{S}(f_j)+\mathit{N}(f_j),j=1,\dots ,J$$

(3)

where $f_L\le f_j\le f_H$ , and $\mathit{A}(f_j,\theta )$ is the $M\times P$ array response matrix

$$\mathit{A}(f_j,\theta )=[\mathit{a}(f_j,{\theta }_1),\mathit{a}(f_j,{\theta }_2),\dots ,\mathit{a}(f_j,{\theta }_P)]$$

(4)

whose columns are the $M\times 1$ steering vectors in which $\mathit{a}(f_j,{\theta }_p)=[1,e^{i2\pi f_{j}d\sin ({\theta }_p)/c},\dots ,e^{i2\pi f_j(M-1)d\sin ({\theta }_p)/c}{]}^T$ is the steering vector at frequency $f_j$ .

For convenience, we use ${\mathit{A}}_j(\theta )$ and ${\mathit{a}}_j({\theta }_p)$ instead of $\mathit{A}(f_j,\theta )$ and $\mathit{a}(f_j,{\theta }_p)$ in the following part of this article. According to equation (3), if the noise is Gaussian white noise with variance ${\sigma }_n^2$ both temporally and spatially, the array covariance matrix at frequency $f_j$ is

$$\mathit{R}(f_j)=E[\mathit{X}(f_j){\mathit{X}}^H(f_j)]={\mathit{A}}_j(\theta ){\mathit{R}}_s(f_j){\mathit{A}}_j^H(\theta )+{\sigma }_n^2\mathit{I}$$

(5)

where ${\mathit{R}}_s(f_j)=E[\mathit{S}(f_j){\mathit{S}}^H(f_j)]$ is the correlation matrix of signals at frequency $f_j$ . If the P signals are uncorrelated, then the $P\times P$ matrix ${\mathit{R}}_s(f_j)$ has full rank. Then, the signal subspace matrix ${\mathit{U}}_s(f_j)$ and the noise subspace matrix ${\mathit{U}}_n(f_j)$ can be obtained from eigenvalue decomposition of the array covariance matrix

$${\mathit{U}}_s(f_j)=[{\mathit{u}}_1(f_j),{\mathit{u}}_2(f_j),\dots ,{\mathit{u}}_P(f_j)]$$

(6)

$${\mathit{U}}_n(f_j)=[{\mathit{u}}_{P+1}(f_j),{\mathit{u}}_{P+2}(f_j),\dots ,{\mathit{u}}_M(f_j)]$$

(7)

where ${\mathit{u}}_i(f_j),i=1,\dots ,M$ are the orthogonal eigenvectors of $\mathit{R}(f_j)$ indexed in descending order according to their corresponding eigenvalues.

Test of orthogonality of subspaces

According to the idea of TOPS in Yoon et al., ⁹ the steering vector ${\mathit{a}}_j({\theta }_j)$ could be transformed into another array steering ${\mathit{a}}_k({\theta }_k)$ through matrix multiplication

$${\mathit{a}}_k({\theta }_k)={\mathit{\Phi }}_q(f_q,{\theta }_q){\mathit{a}}_p({\theta }_p)$$

(8)

where ${\mathit{\Phi }}_q(f_q,{\theta }_q)$ is a diagonal matrix whose elements are $[{\varphi }_q(f_q,{\theta }_q){]}_{m,m}=e^{-i2\pi f_q(m-1)d\sin ({\theta }_q)/c}$ , $f_k=f_p+f_q$ , and $\sin ({\theta }_k)=((f_p/f_k)\sin ({\theta }_p))+((f_p/f_k)\sin ({\theta }_q))$ . If ${\theta }_p={\theta }_q$ , then ${\theta }_k={\theta }_q$ , that is, we can change the steering vector corresponding to one frequency into that of another frequency without changing the DOA.

Define a novel $M\times P$ dimension matrix $\mathit{U}(f_p)$ as follows

$$\mathit{U}(f_p)=\mathit{\Phi }(\Delta f_p,\varphi ){\mathit{U}}_s(f_1),p=2,3,\dots ,J$$

(9)

where $\Delta f_p=f_p-f_1$ , ${\mathit{U}}_s(f_1)$ is the signal subspace at the reference frequency $f_1$ , and $\varphi$ is the possible angle. And then, we can define $P\times (J-1)(M-P)$ matrix $\mathit{D}(\varphi )$

$$\mathit{D}(\varphi )=[{\mathit{U}}^H(f_2){\mathit{U}}_n(f_2),{\mathit{U}}^H(f_3){\mathit{U}}_n(f_3),\dots ,{\mathit{U}}^H(f_J){\mathit{U}}_n(f_J)]$$

(10)

where ${\mathit{U}}_n(f_j)$ is the noise eigenvector associated with the frequency bin $f_j$ . The matrix would lose its rank if $\varphi ={\theta }_p$ for $p=1,2,\dots ,P$ .

TOPS has significant performance in DOA estimation; however, the reference frequency has serious influence on it. When the selected frequency bin is unsuitable, there would appear fake peaks. TOFS ¹¹ constructs the searching steering vectors of every possible DOA and every frequency; thus, it can abstain the fake peaks that are often emerged in TOPS.

For a arbitrarily frequency bin $f_j$ , the following equation is satisfied if $\varphi$ is the true angle of source

$${\mathit{a}}_j^H(\varphi ){\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j){\mathit{a}}_j(\varphi )=0$$

(11)

Then, a new matrix whose element is $U_{j,\varphi }$ can be constructed

$$\mathit{D}\prime (\varphi )=[U_{2,\varphi },\dots ,U_{J,\varphi }]$$

(12)

where $U_{j,\varphi }$ is a scalar which can be denoted as $U_{j,\varphi }={\mathit{a}}_j^H(\varphi ){\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j){\mathit{a}}_j(\varphi ),j=2,3,\dots ,J$ . We can attain the estimated DOA of sources by judging the extent that each element approaches 0, that is, through one-dimensional (1D) angle searching on $\varphi$ , signals’ DOA estimation could be obtained by

$$\hat{\theta }=\arg max\left[\frac{1}{{\sigma }_{min}(\varphi )}\right]$$

(13)

in which ${\sigma }_{min}(\varphi )$ is the smallest singular value of $\mathit{D}\prime (\varphi )$ .

MTOFS

Principle of MTOFS

The TOFS method has good performance; however, it has two intrinsic limitations, one is the heavy computational cost and the other is requiring the number of signals to be known or to be exactly estimated in advance. Generally, the covariance matrix $\mathit{R}(f_j)$ corresponding to $f_j$ is a complex matrix and the inverse of computational complexity is $O(M^3)$ , where M is the number of sensor element. Since multiplication of complex matrix costs four times more than that of real multiplications, and utilizing real matrix instead of complex matrix, we can decrease the whole computational cost of the method.

Suppose we have obtained the prior information of the signals’ number. Then perform eigendecomposition for $\mathit{R}(f_j)$

$$\mathit{R}(f_j)=[\begin{array}{cc}{\mathit{U}}_s(f_j)\hfill & {\mathit{U}}_n(f_j)\hfill \end{array}]\left[\begin{array}{c}\begin{array}{cc}{\mathit{\Sigma }}_s\hfill & \mathbf{0}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\mathbf{0}\hfill & {\sigma }_n^2\mathit{I}\hfill \end{array}}\end{array}\right]\left[\begin{array}{c}{\mathit{U}}_s^H(f_j)\hfill \\ \multicolumn{1}{c}{{\mathit{U}}_n^H(f_j)}\end{array}\right]={\mathit{U}}_s(f_j){\mathit{\Sigma }}_s{\mathit{U}}_s^H(f_j)+{\sigma }_n^2{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)$$

(14)

where ${\mathit{U}}_s(f_j)$ is the signal subspace at frequency $f_j$ , ${\mathit{U}}_n(f_j)$ is the noise subspace, and ${\mathit{\Sigma }}_s=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_P)$ is the diagonal matrix. According to subspace theory, we know

$$\begin{array}{c}{\mathit{U}}_s(f_j){\mathit{U}}_s^H(f_j)+{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)={\mathit{I}}_M\hfill \\ \multicolumn{1}{c}{{\mathit{U}}_s^H(f_j){\mathit{U}}_s(f_j)={\mathit{I}}_P}\\ \multicolumn{1}{c}{{\mathit{U}}_n^H(f_j){\mathit{U}}_n(f_j)={\mathit{I}}_{M-P}}\\ \multicolumn{1}{c}{{\mathit{U}}_s^H(f_j){\mathit{U}}_n(f_j)=\mathbf{0}}\end{array}\}$$

(15)

Performing inverse operation for $\mathit{R}(f_j)$ , we can get

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{R}}^{-1}(f_j)={\mathit{U}}_s(f_j){\mathit{\Sigma }}_s^{-1}{\mathit{U}}_s^H(f_j)+{\sigma }_n^{-2}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)}\\ \multicolumn{1}{c}{={\sigma }_n^{-2}[{\mathit{U}}_s(f_j)\mathrm{diag}\left\{{({\lambda }_i\backslash {\sigma }_n^2)}^{-1}\right\}{\mathit{U}}_s^H(f_j)+{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}\end{array}$$

(16)

where ${\lambda }_i\backslash {\sigma }_n^2$ is the SNR of ith source, when ${\lambda }_i\backslash {\sigma }_n^2>>1$ , $\mathrm{diag}\{({\lambda }_i\backslash {\sigma }_n^2{)}^{-1}\}$ is the zero matrix, and equation (16) can be rewritten as

$${\mathit{R}}^{-1}(f_j)\approx {\sigma }_n^{-2}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)\mathrm{SNR}>>1$$

(17)

Equation (17) testifies that the Capon method has the same performance as the MUSIC method when SNR is large, that is

$$f_{\mathrm{Capon}}(\theta )\approx {\sigma }_n^2{f}_{\mathrm{MUSIC}}(\theta ),\mathrm{SNR}>>1$$

(18)

Since ${\mathrm{Re}}^{-1}(\mathit{R})=2(\mathit{R}+{\mathit{R}}^{*}{)}^{-1}=2{\mathit{R}}^{-1}-2{\mathit{R}}^{-1}(\mathit{I}+{\mathit{R}}^{*}{\mathit{R}}^{-1}{)}^{-1}{\mathit{R}}^{*}{\mathit{R}}^{-1}$ , then

$$\begin{array}{l}{{\mathrm{Re}}^{-1}\left[\mathit{R}(f_j)\right]|}_{\mathrm{SNR}\gg 1}\approx 2{\sigma }_n^{-2}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)-2{\mathit{D}}_1{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)\\ =2({\sigma }_n^{-2}\mathit{I}-{\mathit{D}}_1){\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)\end{array}$$

(19)

where ${\mathit{D}}_1={\mathit{R}}^{-1}(f_j)[\mathit{I}+{\mathit{R}}^{*}(f_j){\mathit{R}}^{-1}(f_j){]}^{-1}{\mathit{R}}^{*}(f_j)$ . Equation (19) indicates that each row of ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ is the row combination of matrix ${\mathit{U}}_n^H(f_j)$ . Since ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ is the Hermitian matrix, ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ is the column combination of matrix ${\mathit{U}}_n(f_j)$ , then

$$\mathrm{span}\{{\mathrm{Re}}^{-1}[\mathit{R}(f_j)]\}\in \mathrm{span}[{\mathit{U}}_n(f_j)],\mathrm{SNR}>>1$$

(20)

Similarly, since ${\mathrm{Re}}^{-1}(\mathit{R})=2({\mathit{R}}^{*}+\mathit{R}{)}^{-1}=2{\mathit{R}}^{-*}-2{\mathit{R}}^{-*}(\mathit{I}+\mathit{R}{\mathit{R}}^{-*}{)}^{-1}\mathit{R}{\mathit{R}}^{-*}$ , then

$$\begin{array}{c}\multicolumn{1}{c}{{{\mathrm{Re}}^{-1}[\mathit{R}(f_j)]|}_{\mathrm{SNR}>>1}\approx 2{\sigma }_n^{-2}{\mathit{U}}_n^{*}(f_j){\mathit{U}}_n^T(f_j)-2{\mathit{D}}_1^{*}{\mathit{U}}_n^{*}(f_j){\mathit{U}}_n^T(f_j)}\\ \multicolumn{1}{c}{=2({\sigma }_n^{-2}\mathit{I}-{\mathit{D}}_1^{*}){\mathit{U}}_n^{*}(f_j){\mathit{U}}_n^T(f_j)}\end{array}$$

(21)

Equation (21) indicates that each row of ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ is also the row combination of matrix ${\mathit{U}}_n^T(f_j)$ . Since ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ is Hermitian matrix, ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ is also the column combination of matrix ${\mathit{U}}_n^{*}(f_j)$ , then

$$\mathrm{span}[R{e}^{-1}\mathit{R}(f_j)]\in \mathrm{span}[{\mathit{U}}_n^{*}(f_j)],\mathrm{SNR}>>1$$

(22)

Combining equations (20) and (22), we know

$$\mathrm{span}\left\{R{e}^{-1}[\mathit{R}(f_j)]\right\}\in \mathrm{span}[{\mathit{U}}_n^{*}(f_j)]\cap \mathrm{span}[{\mathit{U}}_n(f_j)],\mathrm{SNR}>>1$$

(23)

Then, we can define a new matrix

$$\mathit{D}″(\varphi )=[U{\prime }_{2,\varphi },\dots ,U{\prime }_{J,\varphi }]$$

(24)

where $U{\prime }_{j,\varphi }={\mathit{a}}_j^H(\varphi ){\mathrm{Re}}^{-1}[\mathit{R}(f_j)]{\mathit{a}}_j(\varphi ),j=2,3,\dots ,J$ . Note that equation (23) is correct only when the $\mathrm{SNR}>>1$ . However, for the existence of noise, SNR cannot be very large in real world; to exploit equation (23) in practice, we can modify it. We can choose a large appropriate integer coefficient m to weight $\mathit{R}(f_j)$ , that is, using ${\mathit{R}}^m(f_j)$ to replace $\mathit{R}(f_j)$ , then substitute ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ for ${\mathrm{Re}}^{-1}[\mathit{R}(f_j)]$ .

Since

$$\begin{array}{l}{\mathit{R}}^{-m}(f_j)={\left[{\mathit{R}}^{-1}(f_j)\right]}^m={\left[{\mathit{U}}_s(f_j){\Sigma }_s^{-1}{\mathit{U}}_s^H(f_j)+{\sigma }_n^{-2}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)\right]}^m\\ ={\sigma }_n^{-2m}{\mathit{U}}_s(f_j)diag\left\{{({\sigma }_n^2/{\lambda }_i)}^m\right\}{\mathit{U}}_s^H(f_j)+{\sigma }_n^{-2m}{\left[{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)\right]}^m\end{array}$$

(25)

For the value of ${\sigma }_n^2/{\lambda }_i$ is always less than 1, if coefficient m is a larger positive integer, $({\sigma }_n^2\backslash {\lambda }_i{)}^m\approx 0$ will be workable, then $\mathrm{diag}\{({\sigma }_n^2/{\lambda }_i{)}^m\}$ is the zero matrix, and equation (25) can then be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{R}}^{-m}(f_j)={[{\mathit{U}}_s(f_j){\mathit{\Sigma }}_s^{-1}{\mathit{U}}_s^H(f_j)+{\sigma }_n^{-2}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}^m}\\ \multicolumn{1}{c}{\approx {\sigma }_n^{-2m}{[{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}^m}\\ \multicolumn{1}{c}{={\sigma }_n^{-2m}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)}\end{array}$$

(26)

From equation (26), we also know

$$\begin{array}{c}\multicolumn{1}{c}{{[{\mathit{R}}^{-m}(f_j)]}^{*}={\left\{{[{\mathit{U}}_s(f_j){\mathit{\Sigma }}_s^{-1}{\mathit{U}}_s^H(f_j)+{\sigma }_n^{-2}{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}^m\right\}}^{*}}\\ \multicolumn{1}{c}{\approx {\left\{{\sigma }_n^{-2m}{[{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}^m\right\}}^{*}}\\ \multicolumn{1}{c}{={\sigma }_n^{-2m}{[{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j){\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j),\dots ,{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}^{*}}\\ \multicolumn{1}{c}{={\sigma }_n^{-2m}{[{\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)]}^{*}}\\ \multicolumn{1}{c}{={\sigma }_n^{-2m}{\mathit{U}}_n^{*}(f_j){\mathit{U}}_n^T(f_j)}\end{array}$$

(27)

Since

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]={\left\{\frac{{\mathit{R}}^m(f_j)+{[{\mathit{R}}^m(f_j)]}^{*}}{2}\right\}}^{-1}}\\ \multicolumn{1}{c}{=2{\{{\mathit{R}}^m(f_j)+{[{\mathit{R}}^m(f_j)]}^{*}\}}^{-1}}\\ \multicolumn{1}{c}{=2{\mathit{R}}^{-m}(f_j)-2{\mathit{R}}^{-m}(f_j)}\\ \multicolumn{1}{c}{{\{\mathit{I}+{[{\mathit{R}}^m(f_j)]}^{*}{\mathit{R}}^{-m}(f_j)\}}^{-1}{[{\mathit{R}}^m(f_j)]}^{*}{\mathit{R}}^{-m}(f_j)}\\ \multicolumn{1}{c}{=2(\mathit{I}-{\mathit{D}}_2){\mathit{R}}^{-m}(f_j)}\\ \multicolumn{1}{c}{=2{\sigma }_n^{-2m}(\mathit{I}-{\mathit{D}}_2){\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)}\end{array}$$

(28)

where ${\mathit{D}}_2={\mathit{R}}^{-m}(f_j)\{\mathit{I}+[{\mathit{R}}^m(f_j){]}^{*}{\mathit{R}}^{-m}(f_j){\}}^{-1}[{\mathit{R}}^m(f_j){]}^{*}$ ; equation (28) indicates that each row of ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ is the row combination of matrix ${\mathit{U}}_n^H(f_j)$ . Since ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ is the Hermitian matrix, ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ is the column combination of matrix ${\mathit{U}}_n(f_j)$ , then

$$\mathrm{span}\left\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\right\}\in \mathrm{span}[{\mathit{U}}_n(f_j)]$$

(29)

Similarly, since

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]={\left\{\frac{{[{\mathit{R}}^m(f_j)]}^{*}+{\mathit{R}}^m(f_j)}{2}\right\}}^{-1}}\\ \multicolumn{1}{c}{=2{\{{[{\mathit{R}}^m(f_j)]}^{*}+{\mathit{R}}^m(f_j)\}}^{-1}}\\ \multicolumn{1}{c}{=2{[{\mathit{R}}^m(f_j)]}^{-*}-2{[{\mathit{R}}^m(f_j)]}^{-*}}\\ \multicolumn{1}{c}{{\{\mathit{I}+{\mathit{R}}^m(f_j){[{\mathit{R}}^m(f_j)]}^{-*}\}}^{-1}{\mathit{R}}^m(f_j){[{\mathit{R}}^m(f_j)]}^{-*}}\\ \multicolumn{1}{c}{=2(\mathit{I}-{\mathit{D}}_3){[{\mathit{R}}^m(f_j)]}^{-*}}\\ \multicolumn{1}{c}{=2(\mathit{I}-{\mathit{D}}_3){[{\mathit{R}}^m(f_j)]}^{*}}\\ \multicolumn{1}{c}{=2{\sigma }_n^{-2m}(\mathit{I}-{\mathit{D}}_3){\mathit{U}}_n^{*}(f_j){\mathit{U}}_n^T(f_j)}\end{array}$$

(30)

where ${\mathit{D}}_3=[{\mathit{R}}^m(f_j){]}^{-*}\{\mathit{I}+{\mathit{R}}^m(f_j)[{\mathit{R}}^m(f_j){]}^{-*}{\}}^{-1}{\mathit{R}}^m(f_j)$ ; equation (30) indicates that each row of ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ is also the row combination of matrix ${\mathit{U}}_n^T(f_j)$ . Since ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ is the Hermitian matrix, ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ is also the column combination of ${\mathit{U}}_n^{*}(f_j)$ , then

$$\mathrm{span}\left\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\right\}\in \mathrm{span}[{\mathit{U}}_n^{*}(f_j)]$$

(31)

Combining equations (29) and (31), we know

$$\mathrm{span}\left\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\right\}\in \mathrm{span}[{\mathit{U}}_n^{*}(f_j)]\cap \mathrm{span}[{\mathit{U}}_n(f_j)]$$

(32)

Equation (32) indicates that using a large positive integer, $\mathrm{span}\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\}$ can be seen as a subset of $\mathrm{span}[{\mathit{U}}_n(f_j)]$ ; it also contains a part of the vectors of $\mathrm{span}({\mathit{U}}_n(f_j))$ . Therefore, we can utilize $\mathrm{span}\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\}$ instead of $\mathrm{span}[{\mathit{U}}_n(f_j)]$ to estimate signals’ DOAs. Since $\mathrm{span}\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\}$ can approximate to noise subspace very fast, thus overcoming the weakness of equation (23), that is, ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ can substitute ${\mathit{U}}_n(f_j)$ no matter how the SNR is; moreover, no priori information of incident signal number need to be known in the whole process.

Since $\mathrm{span}\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\}$ has a double orthogonality to $\mathrm{span}[{\mathit{A}}_j(\theta )]$ at both the true DOAs and their mirror directions, that is

$$\{\begin{array}{c}\mathrm{span}\left\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\right\}\perp \mathrm{span}[{\mathit{A}}_j(\theta )]\hfill \\ \multicolumn{1}{c}{\mathrm{span}\left\{{\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\right\}\perp \mathrm{span}[{\mathit{A}}_j(-\theta )]}\end{array}$$

(33)

This is because ${\mathit{A}}_j^{*}(\theta )={\mathit{A}}_j(-\theta )$ holds for ULA. Such a merit can help us limit spectral search to only half of the total angular field-of -view.

We now define

$$\mathit{D}‴(\varphi )=[U_{2,\varphi }^{″},\dots ,U_{J,\varphi }^{″}]$$

(34)

where ${U^{\prime }}_{j,\varphi }=\mathrm{Re}[{\mathit{a}}_j^H(\varphi )]\times {\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]\times Re[{\mathit{a}}_j(\varphi )],j=2,3,\dots ,J$ , and m is a larger positive integer.

Via judging the extent that each element approaches 0, we can acquire the DOA estimation of wideband coherent signals by the following equation

$$P_{\mathrm{MTOFS}}(\theta )=\arg max\left[\frac{1}{{\sigma \prime }_{min}(\varphi )}\right]$$

(35)

where $\sigma {\prime }_{min}(\varphi )$ is the minimum singular value of $\mathit{D}‴(\varphi )$ .

Description of the MTOFS

The proposed method can be implemented as follows:

• Step 1: compute the array output data $\mathit{X}(f_j)$ at the frequency bin j via DTFT and compute the covariance matrix $\mathit{R}(f_j)$ by utilizing equation (5).

• Step 2: select suitable coefficient m to weight matrix $\mathit{R}(f_j)$ , and then use ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ to substitute ${\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)$ .

• Step 3: generate $\mathit{D}‴(\varphi )$ for every hypothesized DOA.

• Step 4: construct spectrum function $P_{\mathrm{MTOFS}}(\theta )$ using equation (35), then via peak search from $[0\circ ~90\circ )$ or $(-90\circ ~0\circ ]$ to acquire true DOAs and its mirror. Through following extreme test criterion from $\{{\theta }_p,-{\theta }_p\},p=1,\dots ,P$ , we can remove the mirror sources and acquire true DOAs:

1. If $(1/({\mathit{a}}_j^H({\theta }_p){\mathit{R}}^{-1}(f_j)\mathit{a}({\theta }_p)))>>(1/({\mathit{a}}_j^H(-{\theta }_p){\mathit{R}}^{-1}(f_j)\mathit{a}(-{\theta }_p)))$ , ${\theta }_p$ is true angle;

2. If $(1/({\mathit{a}}_j^H(-{\theta }_p){\mathit{R}}^{-1}(f_j){\mathit{a}}_j(-{\theta }_p)))>>(1/({\mathit{a}}_j^H({\theta }_p){\mathit{R}}^{-1}(f_j){\mathit{a}}_j({\theta }_p)))$ , $-{\theta }_p$ is true angle;

3. If $(1/({\mathit{a}}_j^H({\theta }_p){\mathit{R}}^{-1}(f_j){\mathit{a}}_j({\theta }_p)))\approx (1/({\mathit{a}}_j^H(-{\theta }_p){\mathit{R}}^{-1}(f_j){\mathit{a}}_j^H(-{\theta }_p)))$ , both ${\theta }_p$ and $-{\theta }_p$ are true angles.

Performance analysis of MTOFS

Two main advantages are expected compared with TOFS method, as follows:

1. Based on real-valued calculation, the new method can swiftly obtain the noise subspace, and the computational cost of subspace estimation can be reduced by 75%. Then, using ${\mathrm{Re}}^{-1}[{\mathit{R}}^m(f_j)]$ to substitute ${\mathit{U}}_n(f_j){\mathit{U}}_n^H(f_j)$ , the computational load of the spectrum peak search can be reduce to half of the original; moreover, the calculation of peak search is also based on the real matrix; thus, the whole computational cost is far less than TOFS.

2. Selecting a suitable coefficient m to weight the covariance matrix, the new method can obtain estimated noise subspace without needing the number of signals.

Simulation results

In this section, to verify the performance of the new method, several simulations have been carried out. The array is a ULA with five elements. The sensor spacing is half the wavelength corresponding to the center frequency. The wideband signals are Gaussian processes with zero means and all signals have the same center frequency $f_0=90\mathrm{Hz}$ and the bandwidth $\mathrm{BW}=20\mathrm{Hz}$ . The sampling frequency is $40\mathrm{Hz}$ and the bandwidth is partitioned into $J=11$ frequency bins which can be denoted as $[80,82,84,86,88,90,92,94,96,98,100]\mathrm{Hz}$ ; each frequency bin has 100 snapshots and the Discrete Fourier Transform (DFT) size is 256. The additional noises in all experiments are white Gaussian sources temporally and spatially. The statistical performance is evaluated by performing 100 Monte Carlo runs for each scenario. The average power of signal $p(p=1,2,\dots ,P)$ is $P_p=E[s_p^2(t)]$ . The average power of every signal is the same and the SNR is defined as

$$\mathrm{SNR}=10{\log }_{10}\left(\frac{P_{s_1}}{{\sigma }^2}\right)$$

(36)

in which ${\sigma }^2$ is the variance of the additive Gaussian noise.

The root mean square error (RMSE) of the DOAs’ estimation is defined as

$$\mathrm{RMSE}=\sqrt{\frac{1}{100P}\sum _{i=1}^{100}\sum _{p=1}^P{({\hat{\theta }}_p(i)-{\theta }_p)}^2}$$

(37)

where ${\hat{\theta }}_p(i)$ is the estimate angle for the ith Monte Carlo trial, ${\theta }_p$ is the true angle, and P is the signal number. The probability of resolution is defined as follows: two signals with DOA ${\theta }_1$ and ${\theta }_2$ can be said to be resolved if the estimates ${\hat{\theta }}_1$ and ${\hat{\theta }}_2$ satisfy that both $|{\hat{\theta }}_1-{\theta }_1|$ and $|{\hat{\theta }}_2-{\theta }_2|$ are less than $|{\theta }_1-{\theta }_2|/2$ .

Example 1: 1D spatial spectrum

Two wideband signals are placed at $10\circ$ and $25\circ$ , and the SNR of each signal is 10 dB. We apply the proposed MTOFS to the received signals. For comparison purposes, TOFS are also considered. Figure 1 shows the spatial pseudo-spectrum of TOFS and MTOFS. For each algorithm, the normalization is obtained by dividing the maximum value of its spectrum. The results in Figure 1 indicate that TOFS can precisely estimate the DOAs, that is, it has two sharp spectrum peaks in total angular field-of-view, and MTOFS can also detect signals’ DOAs; however, for the incursion of conjugate noise subspace, it has mirror sources in the symmetrical location of true DOAs. We can utilize this phenomenon to reduce the computational cost of the method just by half search for spectrum, and then just using extreme test which would bring very little computational cost, we can remove the mirror sources and obtain true signals’ DOA. The results of extreme test are listed in Table 1, and since the results of $f({\hat{\theta }}_p)=(1/({\mathit{a}}_j^H({10}^{\circ }){\mathit{R}}^{-1}(f_j){\mathit{a}}_j({10}^{\circ })))$ are far larger than that of $(1/({\mathit{a}}_j^H(-10\circ ){\mathit{R}}^{-1}(f_j){\mathit{a}}_j(-10\circ )))$ , we can draw the conclusion that ${\hat{\theta }}_1=10\circ$ is the true angle and ${\hat{\theta }}_3=-10\circ$ is the mirror; meanwhile, similar conclusion is obtained that ${\hat{\theta }}_2=25\circ$ is the true angle and ${\hat{\theta }}_4=-25\circ$ is the mirror. To verify the efficiency of MTOFS, we use MATLAB 8.0 to compare the simulation time cost by TOFS and MTOFS in Table 2. From Table 2, we can see that the simulation time cost by MTOFS is far less than that of TOFS; therefore, our method is much more efficient than TOFS. Figure 2 demonstrates that the value of coefficient m has significant effect, and MTOFS has better performance if the value of m is larger; however, with the increase in m, the computational burden of the algorithm can be increased sharply; moreover, when m is increased to some extent, the spectrum peak has no obvious change. The results of Figure 2 also confirm that there is no obvious difference when m is 3 or 8; considering efficient implementation of the algorithm, m should be chosen as $3~5$ .

Figure 1.

Space spectrum of TOFS and MTOFS.

Table 1.

Extreme test.

Angle	$10.11\circ$	$25.08\circ$	$-10.13\circ$	$-25.12\circ$
$f({\hat{\theta }}_p)$	245.12	248.12	4.35	4.26

Table 2.

Computer run time ${10}^{-1}\mathrm{s}$ .

TOFS	6.8102	6.7821	6.8211	6.8245
MTOFS	2.4812	2.5324	2.5468	2.5310

TOFS: test of orthogonality of frequency subspaces; MTOFS: modified TOFS.

Figure 2.

Space spectrum of MTOFS with different m.

Example 2: RMSE and probability of resolution versus SNR

In the second case, the RMSE and probability of resolution as a function of SNR are checked. We use ULA with five elements and two wideband signals are impinging from $[10\circ ,20\circ ]$ . The empirical RMSE of MTOFS is shown in Figure 3 and the TOFS is also plotted for comparison. The SNR varies from 0 to 15 dB with 100 snapshots collected in each frequency bin. Figure 3 indicates that the MTOFS has lower RMSE compared with TOFS when SNR is low; however, with the increase in SNR, the RMSE of MTOFS is very close to TOFS and the difference between them is small. To show the influence of the coefficient m for MTOFS, we compare the RMSE of localization under different m with SNR. Figure 4 shows that with the increase in m, MTOFS RMSE is less; however, when m increases to some extent, improvement of MTOFS performance becomes unobvious if we continue to enhance m.

Figure 3.

RMSE versus SNR with TOFS and MTOFS.

Figure 4.

RMSE versus SNR for MTOFS with different m.

We plot the detection probability of resolution curve for the MTOFS and TOFS. The SNR ranges from 0 to 15 dB, and the snapshots of each frequency are fixed at 100. Figure 5 indicates that the detection probability of DOA estimation of the proposed method is smaller than that of TOFS when SNR is lower than 4 dB, while it is almost as good as TOFS when SNR exceeds 4 dB.

Figure 5.

Probability of resolution versus SNR for MTOFS and TOFS.

Example 3: RMSE and probability of resolution versus snapshots

To further check the performance of the MTOFS, the RMSE of the estimated parameters versus snapshots are shown in Figure 6. The number of snapshots varies from 50 to 1000, the SNR is fixed at 10 dB, and the other conditions are similar to those in Example 2.

Figure 6.

RMSE versus snapshots with MTOFS and TOFS.

The empirical RMSE of MTOFS is shown in Figure 6, and it is compared with TOFS. We can see that MTOFS is superior to the TOFS when the snapshots are small due to the fact that the rank loss of complex covariance matrix is more serious than the real covariance matrix, and with the increase in the snapshots, the difference between them is almost indistinguishable. Figure 7 confirms the effectiveness of the coefficient m, and the results are displayed in Figure 7, from which we can see that the results are similar to those in Figure 4, that is, MTOFS RMSE is less when m is larger. Figure 8 indicates that the detection probability of DOA estimation of the proposed method is a little bit inferior to the TOFS; however, with the augment of snapshots, the difference between them can be neglected.

Figure 7.

RMSE versus snapshots for MTOFS with different m.

Figure 8.

Probability of resolution versus snapshots for MTOFS and TOFS.

Conclusion

In this article, we propose an efficient DOA estimation method for wideband signals based on ULA. The first favorable merit of the proposed algorithm is that it does not need to know the source number information, and such a merit is extremely desirable in real-world applications since the detection of source number is a highly tough task. The second merit is that the proposed method is a low computational complexity algorithm and the computational cost is far less than that of TOFS. The simulation results demonstrate that the proposed method is more valuable and effective than TOFS.