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

Authors: Shuang Li; Xiaoxiao Jiang; Sai Ma; Yingguan Wang
Journal: International Journal of Antennas and Propagation
 (2013)

---

## 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.


---

## Body

## 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 ℓ0 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ℓ1 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 2M-2 sources by using a uniform linear array (ULA) with M elements. Moreover, our method does not have to know the number of sources in advance.The paper is organized as follows. Section2 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(·)*, (·)T, (·)H, and E{·} denote complex conjugation, matrix transposition, matrix complex conjugate plus transposition, and the statistical expectation, respectively.Consider this scenario thatQ narrowband far-field sources impinge on a ULA with M sensors corrupted by additive Gaussian noise. The observed signal at the mth sensor is
(1)xm(t)=∑i=1Qamisi(t)+nm(t),m=1,2,…,M,
where
(2)ami=exp(j2πλdmsinθi),dm=(m-1)d denotes the distance from the mth sensor to the reference sensor, d denotes the intersensor spacing, λ denotes the carrier wavelength, si(t) denotes the ith source signal, θi denotes the DOA of the ith signal, and nm(t) denotes the received noise of the mth sensor.Letx(t)=[x1(t),x2(t),…,xM(t)]T denote the received signal vector, A(θ) denote the array manifold matrix, whose (m,i)th entry is ami, s(t)=[s1(t),s2(t),…,sQ(t)]T denote the signal vector, and n(t)=[n1(t),n2(t),…,nM(t)]Tdenote the noise vector. Then, (1) can be expressed in a compact form as
(3)x(t)=A(θ)s(t)+n(t).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]:
(4)cum(xk(t),xl*(t),xp(t),xq*(t))=E{xk(t)xl*(t)xp(t)xq*(t)}-E{xk(t)xl*(t)}E{xp(t)xq*(t)}-E{xk(t)xp(t)}E{xl*(t)xq*(t)}-E{xk(t)xq*(t)}E{xl*(t)xp(t)},hhhhhhhhhhhhhhhl1≤k,l,p,q≤M.
Under the assumptions (A2) and (A3) and by using the properties of cumulants,
(5)cum(xk(t),xl*(t),xp(t),xq*(t))=∑i=1Qakiali*apiaqi*cum(si,si*,si,si*)hhhhhhhhhhhhh1≤k,l,p,q≤M.Note that the noise term is vanished since the fourth order cumulant of Gaussian noise is zero.Substituting (2) into (5), we have
(6)cum(xk(t),xl*(t),xp(t),xq*(t))=∑i=1Qexp(j2πλ(dk-dl+dp-dq)sin(θi))×cum(si,si*,si,si*)=∑i=1Qexp(j2πλ(k+p-l-q)dsin(θi))×cum(si,si*,si,si*)=b~k+p-l-qTcs,1≤k,l,p,q≤M,
where b~k+p-l-q=[exp(j(2π/λ)(k+p-l-q)dsinθ1),…,exp(j(2π/λ)(k+p-l-q)dsinθQ)]T and cs=[cum(s1,s1*,s1,s1*),cum(s2,s2*,s2,s2*),…,cum(sQ,sQ*,sQ,sQ*)]T. When k, l, p, and q vary from 1 to M, k+p-l-q takes values from 2(1-M) to 2(M-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 k+p-l-q, giving rise to a column vector c.  Define a matrix B∈ℂ(4M-3)×Q whose (m,i)th entry is exp(j(2π/λ)(m-(2M-1))dsinθi). Then we have
(7)c=Bcs.
Comparing (7) to (3), we can see that c behaves like the received signal vector at an array whose array manifold matrix is given by B and the source signal vector is given by cs. When the sensor location of the physical array is given at {(k-1)d,1≤k≤M}, the sensor location of the virtual array is {(k+p-l-q)d,1≤k,p,l,q≤M}.We now sparsely represent the new signal vectorc in a redundant basis. Define a set Π={θ-1,θ-2,…,θ-Nθ}, which denotes potential source locations of interest. The number of the potential source locations Nθ should be much greater than the number of actual sources Q and the number of sensors M. Define the overcomplete basis B-=[b(θ-1),b(θ-2),…,b(θ-Nθ)] and the spatial spectrum ps=[P1,P2,…,PNθ], where b(θ) denotes the steering vector of the virtual array and Pj=cs(i) if θ-j=θi; otherwise Pj=0. Thus we obtain
(8)c=B-ps.
In practice, the vector c is estimated by (4). We assume
(9)c^=c+Δc,
where c^ denotes the estimate value of c and Δc denotes the estimate errors. Substituting (9) into (8), we have
(10)c^=B-ps+Δc.
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:
(11)minps∥ps∥1,subjectto∥c^-B-ps∥22≤β.β is also termed regularization parameter, which trades off the sparse fitting noise (the ℓ2 norm) and the sparsity of ps. 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 ps, 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 vectorc^. 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 O(TM4) and solving (11) requires O(Nθ3), where T 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 O(Q3Nθ3), which grows proportionally with the number of sources Q. Obviously, the computational cost of NR-SpCF is somewhat higher than L1-SVD and L1-SRACV if Q=1. However since Nθ is much greater than M, the computational cost is lower than L1-SVD and L1-SRACV when Q>1.  For example, T=1000, M=6, Q=3, and Nθ=180; then TM4=1.296×106, Nθ3≈5.8×106, Q3Nθ3≈1.57×108. 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 O(Nθ3). 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ΦX=Z, where Φ, X, and Z denote the measured matrix, the sparse matrix, and the observed matrix, respectively. With the assumption that any h columns of Φ are linearly independent and rank(Z)=L<h, a solution with number of nonzero entries r, where r≤ru=⌈(h+L)/2⌉-1, is unique (where ⌈·⌉ denotes the ceiling operation) [29].In NR-SpCF,rank(c^)=L=1, h=4M-3, the maximal number of sources that can be resolved is ru=⌈(4M-3+1)/2⌉-1=2M-2. It is equivalent to that of 4-MUSIC for a ULA with M 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(xk(t),xl*(t), xp(t),xq*(t))} in ascending order of k+p-l-q, forming the vector c^.(S3)
Solve (11) by using the software package CVX or SeDuMi to obtain the sparse spatial spectra ps.(S4)
Find the maxima ofps, 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 M=6, Q=2, θ={-6.2°,25.6°}, SNR=10 dB, and T=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.Spatial spectra of NR-SpCF for differentβ.
(a)
(b)
(c)In last section, a lemma was introduced to interpret that the proposed method NR-SpCF can resolve up to2M-2 sources using a ULA with M sensors. Now we give an experiment to verify it. In this example, we set M=4, Q=6, T=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, whenM=4, Q=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(12)RMSE=1QNmc∑q=1Q∑i=1Nmc(θ^q,i-θq,i),
where Nmc denotes the number of Monte Carlo trials and θ^q,i and θq,i denote the estimate value and the real value of the qth signal in the ith trial. Figure 3 illustrates the RMSE as a function of SNR for the six methods by averaging 100 independent trials, where M=6, Q=2, θ={4.5°,25.8°}, and T=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.M=6, Q=2, and T=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 Figure4, 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.M=6, Q=2, and T=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 beQ=2, M=6, θ={4.5°,25.8°}, SNR=5 dB, and Nmc=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 |θ^1-θ1| and |θ^2-θ2| are smaller than |θ1-θ2|/2, where θ1 and θ2 denote the true DOAs and θ^1 and θ^2 denote the estimating DOAs. The probabilities of resolution are illustrated in Figure 6, where M=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.M=6, and T=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.



---

*Source: 897158-2013-12-16.xml*