Sparse Cumulants Fitting for Direction-of-Arrival Estimation without Redundancy

Abstract

A novel direction-of-arrival (DOA) estimation method is proposed based on the sparse cumulants fitting without redundancy. Firstly, we derive that some fourth order cumulants of the array output are redundant and therefore are removed to reduce computational complexity. Then, the left cumulants are sparsely represented on an overcomplete basis and the DOAs are resolved by using a software package. Despite introducing a high variance, the proposed method shows several advantages including the ability to detect more sources than sensors, high resolution, and robustness to all kinds of Gaussian noise. Besides, our method does not have to know, a priori, the number of sources. Simulation results are presented to illustrate the effectiveness and efficiency of the proposed method.

1. Introduction

Direction-of-arrival (DOA) estimation using sensor arrays is an important topic since it has been widely used in many applications such as radar, sonar, and wireless communications. Plenty of methods [1–10] were proposed in decades. Among these, the subspace-based methods such as MUSIC [2] and ESPRIT [3] were most popular, which show superior performance than beamforming [4] and time-difference-estimation-based methods [5]. By utilizing intelligent algorithms such as support vector regression and neural network, excellent performance was also achieved for DOA estimation [6–8].

Recently, thanks to the development of compressive sensing (CS) theories [11, 12], the DOA estimation problem has attracted considerable attention. Malioutov et al. [13] proposed the L1-SVD method, which is one of the most successful CS-based DOA estimation methods, showing superior performance than conventional methods. Hyder and Mahata [14] minimized the mixed norm approximation and proposed a method called Joint ℓ₀ approximation (JZLA). By vectorizing the covariance matrix of the array output, Zeng et al. [15–17] presented a sparse spectral fitting (SpSF) method which can estimate the DOAs as well as the power at each DOA. Stocia et al. [18] proposed an iterative method, termed SPICE, which is based on covariance matrix fitting and avoids selecting regularization parameter. By using a sparse representation of the array covariance vectors, Yin and Chen [19] presented a method called L1-SRACV and a new error-suppression criterion was given based on weighted covariance vectors fitting.

However, all the methods mentioned above either exploit the sources directly or are based on the second order statistics. In fact, most of the source signals are non-Gaussian, so higher order cumulants (HOC) of the signals contain much more useful information. More importantly, the performance of these methods degrades significantly in spatially correlated noise case. It has been shown that by employing the HOC one cannot only suppress any kind of Gaussian noise but also detect more sources than sensors for DOA estimation [20–23].

Another disadvantage of the above methods is that they have to determine the number of sources before estimating the DOAs. The most commonly used way to determine it is to make use of the information theoretic criteria, like the minimum description length (MDL) or Akaike information criteria (AIC) [24]. However, the probability of successful detection is rather low even in moderate signal-to-noise ratio (SNR) [25]. The failure to determine the number of sources significantly degrades the performance of the above DOA estimation methods.

In this paper, we propose a novel DOA estimation method based on sparse representation of fourth order cumulants. The cumulants are not exploited directly; instead, they are sparsely represented through removing the redundant items from some of cumulants. By enforcing sparsity using ℓ₁ norm just like the method L1-SVD, the spatial spectrum is given by finding the sparsest solution in a redundant basis. For short, the proposed method is abbreviated as NR-SpCF. Despite introducing a high variance, NR-SpCF shows high resolution and robustness to Gaussian noise and can resolve 2𝑀−2 sources by using a uniform linear array (ULA) with 𝑀 elements. Moreover, our method does not have to know the number of sources in advance.

The paper is organized as follows. Section 2 introduces the data model and some assumptions. The details of the proposed method are presented in Section 3. It is followed by simulation results in Section 4. At last, conclusions are drawn in Section 5.

2. Data Model and Assumptions

First, we declare that, in the remainder of the paper, the symbols (⋅)^∗, (⋅)^𝑇, (⋅)^𝐻, and 𝐸{⋅} denote complex conjugation, matrix transposition, matrix complex conjugate plus transposition, and the statistical expectation, respectively.

Consider this scenario that 𝑄 narrowband far-field sources impinge on a ULA with 𝑀 sensors corrupted by additive Gaussian noise. The observed signal at the 𝑚th sensor is 𝑥_𝑚(𝑡)=^𝑄_𝑖=1𝑎_𝑚𝑖𝑠_𝑖(𝑡)+𝑛_𝑚(𝑡),𝑚=1,2,…,𝑀, (1) where 𝑎_𝑚𝑖=exp𝑗2𝜋𝜆𝑑_𝑚sin𝜃_𝑖, (2)𝑑_𝑚=(𝑚−1)𝑑 denotes the distance from the 𝑚th sensor to the reference sensor, 𝑑 denotes the intersensor spacing, 𝜆 denotes the carrier wavelength, 𝑠_𝑖(𝑡) denotes the 𝑖th source signal, 𝜃_𝑖 denotes the DOA of the 𝑖th signal, and 𝑛_𝑚(𝑡) denotes the received noise of the 𝑚th sensor.

Let 𝐱(𝑡)=[𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑀(𝑡)]^𝑇 denote the received signal vector, 𝐀(𝜽) denote the array manifold matrix, whose (𝑚,𝑖)th entry is 𝑎_𝑚𝑖, 𝐬(𝑡)=[𝑠₁(𝑡),𝑠₂(𝑡),…,𝑠_𝑄(𝑡)]^𝑇 denote the signal vector, and 𝐧(𝑡)=[𝑛₁(𝑡),𝑛₂(𝑡),…,𝑛_𝑀(𝑡)]^𝑇denote the noise vector. Then, (1) can be expressed in a compact form as 𝐱(𝑡)=𝐀(𝜽)𝐬(𝑡)+𝐧(𝑡). (3)

Regarding the sources and the noise, we make the following assumptions.(A1)The sources are non-Gaussian.(A2)The sources are independent of each other and not correlated with the noise.(A3)The noise is Gaussian spatially and temporally white or colored.

3. NR-SpCF: Nonredundant Sparse Cumulants Fitting Method

The definition of fourth order cumulant is as follows [25]: cum𝑥_𝑘(𝑡),𝑥^∗_𝑙(𝑡),𝑥_𝑝(𝑡),𝑥^∗_𝑞(𝑡𝑥)=𝐸_𝑘(𝑡)𝑥^∗_𝑙(𝑡)𝑥_𝑝(𝑡)𝑥^∗_𝑞(𝑡𝑥)−𝐸_𝑘(𝑡)𝑥^∗_𝑙(𝑡)𝐸𝑥_𝑝(𝑡)𝑥^∗_𝑞(𝑡𝑥)−𝐸_𝑘(𝑡)𝑥_𝑝(𝑡)𝐸𝑥^∗_𝑙(𝑡)𝑥^∗_𝑞(𝑡𝑥)−𝐸_𝑘(𝑡)𝑥^∗_𝑞(𝑡)𝐸𝑥^∗_𝑙(𝑡)𝑥_𝑝(𝑡),1≤𝑘,𝑙,𝑝,𝑞≤𝑀. (4) Under the assumptions (A2) and (A3) and by using the properties of cumulants, cum𝑥_𝑘(𝑡),𝑥^∗_𝑙(𝑡),𝑥_𝑝(𝑡),𝑥^∗_𝑞(𝑡=_𝑄𝑖=1)𝑎_𝑘𝑖𝑎^∗_𝑙𝑖𝑎_𝑝𝑖𝑎^∗_𝑞𝑖cum𝑠_𝑖,𝑠^∗_𝑖,𝑠_𝑖,𝑠^∗_𝑖1≤𝑘,𝑙,𝑝,𝑞≤𝑀. (5)

Note that the noise term is vanished since the fourth order cumulant of Gaussian noise is zero.

Substituting (2) into (5), we have cum𝑥_𝑘(𝑡),𝑥^∗_𝑙(𝑡),𝑥_𝑝(𝑡),𝑥^∗_𝑞(𝑡=_𝑄𝑖=1)exp𝑗2𝜋𝜆𝑑_𝑘−𝑑_𝑙+𝑑_𝑝−𝑑_𝑞sin𝜃_𝑖×cum𝑠_𝑖,𝑠^∗_𝑖,𝑠_𝑖,𝑠^∗_𝑖=^𝑄_𝑖=1exp𝑗2𝜋𝜆𝑘+𝑝−𝑙−𝑞𝑑sin𝜃_𝑖×cum𝑠_𝑖,𝑠^∗_𝑖,𝑠_𝑖,𝑠^∗_𝑖=𝐛^𝑇_{𝑘+𝑝−𝑙−𝑞}𝐜_𝑠,1≤𝑘,𝑙,𝑝,𝑞≤𝑀, (6) where 𝐛_{𝑘+𝑝−𝑙−𝑞}=[exp(𝑗(2𝜋/𝜆)(𝑘+𝑝−𝑙−𝑞)𝑑sin𝜃₁),…,exp(𝑗(2𝜋/𝜆)(𝑘+𝑝−𝑙−𝑞)𝑑sin𝜃_𝑄)]^𝑇 and 𝐜_𝑠=[cum(𝑠₁,𝑠^∗₁,𝑠₁,𝑠^∗₁),cum(𝑠₂,𝑠^∗₂,𝑠₂,𝑠^∗₂),…,cum(𝑠_𝑄,𝑠^∗_𝑄,𝑠_𝑄,𝑠^∗_𝑄)]^𝑇. When 𝑘, 𝑙, 𝑝, and 𝑞 vary from 1 to 𝑀, 𝑘+𝑝−𝑙−𝑞 takes values from 2(1−𝑀) to 2(𝑀−1) and some values repeat. Therefore some cumulants are equivalent and redundant. Once these redundant cumulants are removed, the computational complexity is naturally reduced. Then the left items are sorted with ascending order of 𝑘+𝑝−𝑙−𝑞, giving rise to a column vector 𝐜. Define a matrix 𝐁∈ℂ^{(4𝑀−3)×𝑄} whose (𝑚,𝑖)th entry is exp(𝑗(2𝜋/𝜆)(𝑚−(2𝑀−1))𝑑sin𝜃_𝑖). Then we have 𝐜=𝐁𝐜_𝑠. (7) Comparing (7) to (3), we can see that 𝐜 behaves like the received signal vector at an array whose array manifold matrix is given by 𝐁 and the source signal vector is given by 𝐜_𝑠. When the sensor location of the physical array is given at {(𝑘−1)𝑑,1≤𝑘≤𝑀}, the sensor location of the virtual array is {(𝑘+𝑝−𝑙−𝑞)𝑑,1≤𝑘,𝑝,𝑙,𝑞≤𝑀}.

We now sparsely represent the new signal vector 𝐜 in a redundant basis. Define a set Π={𝜃₁,𝜃₂,…,𝜃_𝑁𝜃}, which denotes potential source locations of interest. The number of the potential source locations 𝑁_𝜃 should be much greater than the number of actual sources 𝑄 and the number of sensors 𝑀. Define the overcomplete basis 𝐁=[𝐛(𝜃₁),𝐛(𝜃₂),…,𝐛(𝜃_𝑁𝜃)] and the spatial spectrum 𝐩_𝑠=[𝑃₁,𝑃₂,…,𝑃_𝑁𝜃], where 𝐛(𝜃) denotes the steering vector of the virtual array and 𝑃_𝑗=𝐜_𝑠(𝑖) if 𝜃_𝑗=𝜃_𝑖; otherwise 𝑃_𝑗=0. Thus we obtain 𝐜=𝐁𝐩_𝑠. (8) In practice, the vector 𝐜 is estimated by (4). We assume 𝐜̂=𝐜+Δ𝐜, (9) where 𝐜̂ denotes the estimate value of 𝐜 and Δ𝐜 denotes the estimate errors. Substituting (9) into (8), we have 𝐜̂=𝐁𝐩_𝑠+Δ𝐜. (10) To model the estimate noise well, we introduce a parameter 𝛽. Under the assumption of sparse spatial spectrum, the DOA estimation problem is formulated to the following convex optimization problem: min_𝐩𝑠𝐩_𝑠1,subjectto𝐜̂−𝐁𝐩_𝑠²₂≤𝛽. (11)𝛽 is also termed regularization parameter, which trades off the sparse fitting noise (the ℓ₂ norm) and the sparsity of 𝐩_𝑠. Generally, 𝛽 is empirically selected. The optimization problem can be resolved by using some software packages CVX [26] or SeDuMi [27]. After finding the spatial spectrum 𝐩_𝑠, the DOAs can be obtained by finding the maxima of the spatial spectrum.

Remark 1. Since 𝛽 has great impact on the solution of (11), it is of utmost importance to choose the right 𝛽. Although under a few assumptions some methods were given in [19, 28], generally there is no proper way to choose it until now. And it is still an open problem about how to select the regularization parameter. An example will be given in the simulation section to show the importance of selecting a proper 𝛽. In all the simulation tests, the regularization parameter will be selected in advance manually.

Remark 2. All the cumulants are not used to construct the vector 𝐜̂. We only choose one element if there are multiplicities at a given virtual sensor location. Some cumulants corresponding to the redundant sensors are removed to reduce the computational complexity. Regarding the computational cost of NR-SpCF, calculating of the fourth order cumulants requires 𝑂(𝑇𝑀⁴) and solving (11) requires 𝑂(𝑁³_𝜃), where 𝑇 denotes the number of snapshots. The computational cost of L1-SVD [13] and L1-SRACV [19] is mainly concentrated on solving an optimization problem and is 𝑂(𝑄³𝑁³_𝜃), which grows proportionally with the number of sources 𝑄. Obviously, the computational cost of NR-SpCF is somewhat higher than L1-SVD and L1-SRACV if 𝑄=1. However since 𝑁_𝜃 is much greater than 𝑀, the computational cost is lower than L1-SVD and L1-SRACV when 𝑄>1. For example, 𝑇=1000, 𝑀=6, 𝑄=3, and 𝑁_𝜃=180; then 𝑇𝑀⁴=1.296×10⁶, 𝑁³_𝜃≈5.8×10⁶, 𝑄₃𝑁³_𝜃≈1.57×10⁸. It can be clearly seen that the computational cost of NR-SpCF is lower than that of L1-SVD and L1-SRACV. However, it is higher than that of 4-MUSIC [22], where the main cost is in calculation of fourth order cumulants and eigen decomposition of a matrix, and SpSF [15], where the computational cost is mainly 𝑂(𝑁³_𝜃). Although the computational complexity of NR-SpCF is somewhat higher than 4-MUSIC and SpSF, the advantages that we obtain include ability to detect more sources than sensors, high resolution and robustness to all kinds of Gaussian noise.

It is an interesting question of how many sources NR-SpCF can detect. Here we introduce a lemma to determine the number of sources which can be resolved.

Lemma 3. Consider a multiple measurement vectors problem of Φ𝐗=𝐙, where Φ, 𝐗, and 𝐙 denote the measured matrix, the sparse matrix, and the observed matrix, respectively. With the assumption that any ℎ columns of Φ are linearly independent and rank(𝐙)=𝐿<ℎ, a solution with number of nonzero entries 𝑟, where 𝑟≤𝑟_𝑢=⌈(ℎ+𝐿)/2⌉−1, is unique (where ⌈⋅⌉ denotes the ceiling operation) [29].

In NR-SpCF, rank(𝐜̂)=𝐿=1, ℎ=4𝑀−3, the maximal number of sources that can be resolved is 𝑟_𝑢=⌈(4𝑀−3+1)/2⌉−1=2𝑀−2. It is equivalent to that of 4-MUSIC for a ULA with 𝑀 elements. In addition, by using some particular arrays such as minimum-redundancy linear arrays [30] and nested array [31], NR-SpCF can resolve more sources.

Now we summarize our method as follows.(S1)Estimate all the cumulants using (4).(S2)Remove the redundant cumulants and arrange the left cumulants {cum(𝑥_𝑘(𝑡),𝑥^∗_𝑙(𝑡), 𝑥_𝑝(𝑡),𝑥^∗_𝑞(𝑡))} in ascending order of 𝑘+𝑝−𝑙−𝑞, forming the vector 𝐜̂.(S3)Solve (11) by using the software package CVX or SeDuMi to obtain the sparse spatial spectra 𝐩_𝑠.(S4)Find the maxima of 𝐩_𝑠, and the DOAs are obtained.

4. Simulation Results

In this section, some numerical tests are presented to illustrate the performance of NR-SpCF. It is compared with several existing methods including L1-SVD [13], L1-SRACV [19], SpSF [15], MUSIC [2], and 4-MUSIC [22]. We consider a case in which several sources impinge on a ULA separated by half of the wavelength. The sources are modeled as far-field narrowband exponential process. The noise is assumed to be Gaussian white noise unless otherwise noted. The sample grid of potential sources is uniform with 0.1 degrees. We assume that all the methods know the number of sources in advance except NR-SpCF.

First we give an example to illustrate how important it is to select the right regularization parameter 𝛽. Here the parameters are set to be 𝑀=6, 𝑄=2, 𝜃={−6.2^∘,25.6^∘}, SNR=10 dB, and 𝑇=1000. The spatial spectra of NR-SpCF are depicted in Figure 1 for three different values of 𝛽. It can be seen that spurious peaks appear when 𝛽 is too small and some peaks disappear when 𝛽 is too large. The choice of regularization parameter 𝛽 without the information of the noise is still an open problem [13]. Since it is of great importance to select 𝛽, in all the following examples, 𝛽 is selected manually in advance.

(a)

(b)

(c)

Figure 1 

Spatial spectra of NR-SpCF for different 𝛽.

In last section, a lemma was introduced to interpret that the proposed method NR-SpCF can resolve up to 2𝑀−2 sources using a ULA with 𝑀 sensors. Now we give an experiment to verify it. In this example, we set 𝑀=4, 𝑄=6, 𝑇=3000, and 𝜽={−50^∘,−33^∘,−16^∘,1^∘,18^∘,35^∘}. Since L1-SVD, L1-SRACV, SpSF, and MUSIC cannot work in underdetermined case, we compare NR-SpCF to 4-MUSIC and depict the spatial spectra of the two methods in Figure 2. It can be noted that both of the two methods can resolve the 6 sources well and NR-SpCF shows sharper peaks compared with 4-MUSIC. It can be deduced that the resolution of NR-SpCF outperforms that of 4-MUSIC, which will be investigated later.

Figure 2 

Spatial spectra of NR-SpCF and 4-MUSIC for underdetermined DOA estimation case, when 𝑀=4, 𝑄=6. The real DOAs are marked by the circles.

We now evaluate the root square mean error (RMSE) of NR-SpCF as a function of SNR through Monte Carlo simulation. The RMSE is defined as RMSE=1𝑄𝑁_mc^𝑄𝑁_𝑞=1^mc_𝑖=1𝜃_𝑞,𝑖−𝜃_𝑞,𝑖, (12) where 𝑁_mc denotes the number of Monte Carlo trials and 𝜃_𝑞,𝑖 and 𝜃_𝑞,𝑖 denote the estimate value and the real value of the 𝑞th signal in the 𝑖th trial. Figure 3 illustrates the RMSE as a function of SNR for the six methods by averaging 100 independent trials, where 𝑀=6, 𝑄=2, 𝜽={4.5^∘,25.8^∘}, and 𝑇=1000. As we can see, among these methods, L1-SRACV and MUSIC have the best performance, while NR-SpCF shows a little higher RMSE compared to other methods. It is most probably in that the large estimate error of cumulants results in a relative high RMSE of NR-SpCF. Despite its large errors, NR-SpCF obtains many merits including high resolution, capability to detect more sources than sensors, robustness to spatially correlated noise, and Nonnecessity to determine the number of sources in advance.

Figure 3 

RMSE versus SNR for the six methods in Gaussian white noise case. 𝑀=6, 𝑄=2, and 𝑇=1000.

In order to verify that NR-SpCF is insensitive to Gaussian spatially correlated noise, we depict the RMSE of the six methods versus SNR for spatially correlated noise in Figure 4, where the parameters are kept the same as before. The directional-noise is generated using the method in [20] and is uniformly distributed at the area [−30^∘,30^∘]. From Figure 4, we can see that, compared with other methods, NR-SpCF and 4-MUSIC perform a little worse in high SNR and, however, show better performance in low SNR. The former is the reason that, in high SNR, the noise has less impact on these methods, and the latter is because NR-SpCF and 4-MUSIC are based on fourth order cumulants and are blind to Gaussian spatially correlated noise.

Figure 4 

RMSE versus SNR for the six methods in spatially correlated noise case. 𝑀=6, 𝑄=2, and 𝑇=1000.

Since the estimates of cumulants are not accurate and strongly correlated to the number of snapshots, it is necessary to investigate the RMSE of NR-SpCF against the number of snapshots. In this experiment, we set the parameters to be 𝑄=2, 𝑀=6, 𝜽={4.5^∘,25.8^∘}, SNR=5 dB, and 𝑁_mc=100. The RMSE of the six methods as a function of the number of snapshots is given in Figure 5. The results show that NR-SpCF achieves a desirable performance when the number of snapshots is greater than 500. Despite a relative large error of estimating the cumulants, the optimization problem shows somewhat noise-suppression capability, which makes NR-SpCF achieve a desirable performance at a relative small number of samples. In addition, the number of snapshots has a great impact on L1-SRACV. When the number of snapshots is fewer than 200, the RMSE is significantly increased since L1-SRACV is based on the covariance vectors of the array output and the weight the method used is also correlated to the covariance matrix. On the contrary, the RMSE of L1-SRACV is close to MUSIC when the number of snapshots is large enough in that L1-SRACV introduces an error-suppression rule. L1-SVD shows a little change when the number of snapshots is different since this method exploits the observed signals directly.

Figure 5 

The RMSE versus the number of snapshots for the six methods.

In this test, the probabilities of resolution are calculated. We consider two closely spaced signals at {3.8^∘,6.8^∘}. By the definition in [32], the two sources are identified in a trial if both |𝜃₁−𝜃₁| and |𝜃₂−𝜃₂| are smaller than |𝜃₁−𝜃₂|/2, where 𝜃₁ and 𝜃₂ denote the true DOAs and 𝜃₁ and 𝜃₂ denote the estimating DOAs. The probabilities of resolution are illustrated in Figure 6, where 𝑀=6, SNR=5 dB, and 100 Monte Carlo trials are carried out. It can be noted that NR-SpCF achieves higher resolution than the other methods. Both the extended array aperture and the use of sparse signal recovery method make NR-SpCF achieve higher resolution than other methods.

Figure 6 

The resolution ability versus SNR for the six methods. 𝑀=6, and 𝑇=1000.

5. Conclusions

In this paper, a novel DOA estimation method is proposed based on sparse cumulants fitting. We remove the redundant cumulants and sparsely represent the left cumulants in an overcomplete basis. Then the DOAs are resolved by using a software package. In spite of a little higher RMSE, our method achieves higher resolution and can identify more sources than sensor. Moreover, our method shows better performance for spatially correlated noise case. Besides, our method does not require knowing the number of signals in advance. Although our method is based on fourth order cumulants, our method shows moderate complexity. Future research includes wideband DOA estimation based on sparse representation of four order cumulants.

Conflict of Interests

There is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported in part by the Strategic Priority Research Program of the Chinese Academy of Sciences under Grant no. XDA06020300.