Joint Smoothed l₀-Norm DOA Estimation Algorithm for Multiple Measurement Vectors in MIMO Radar

Abstract

Direction-of-arrival (DOA) estimation is usually confronted with a multiple measurement vector (MMV) case. In this paper, a novel fast sparse DOA estimation algorithm, named the joint smoothed $l_0$-norm algorithm, is proposed for multiple measurement vectors in multiple-input multiple-output (MIMO) radar. To eliminate the white or colored Gaussian noises, the new method first obtains a low-complexity high-order cumulants based data matrix. Then, the proposed algorithm designs a joint smoothed function tailored for the MMV case, based on which joint smoothed $l_0$-norm sparse representation framework is constructed. Finally, for the MMV-based joint smoothed function, the corresponding gradient-based sparse signal reconstruction is designed, thus the DOA estimation can be achieved. The proposed method is a fast sparse representation algorithm, which can solve the MMV problem and perform well for both white and colored Gaussian noises. The proposed joint algorithm is about two orders of magnitude faster than the $l_1$-norm minimization based methods, such as $l_1$-SVD (singular value decomposition), RV (real-valued) $l_1$-SVD and RV $l_1$-SRACV (sparse representation array covariance vectors), and achieves better DOA estimation performance.

1. Introduction

Colocated multiple-input multiple-output (MIMO) radar has attracted a growing interest recently because it can achieve higher resolution and better parameter identification compared with conventional phased-array radar [1]. As one type of colocated MIMO radar, the monostatic MIMO radar is equipped with closely-located transmit and receive arrays, which result in the same angle for direction-of-departure (DOD) and direction-of-arrival (DOA) [2]. Aiming at DOA estimation, a large quantity of methods have been proposed, most of them are based on the signal and noise subspaces [3,4,5,6], such as the representative MUSIC (multiple signal classification) algorithm [3] and ESPRIT (estimation of signal parameters via rotational invariance techniques) algorithm [4].

In the area of sensor array signal processing, sparse representation is an important technique to estimate the parameters. When reconstructing the sparse signal, to avoid the NP-hard (non-deterministic polynomial-time hard) $l_0$-norm minimization, different methods such as those based on the relaxed constraint $l_1$-norm minimization [7], the focal underdetermined system solution (FOCUSS) [8] and the sparse Bayesian learning (SBL) [9] have been proposed. These algorithms were originally designed for the single measurement vector (SMV) problem. Then, to extend them to the multiple measurement vector (MMV) case, some algorithms, e.g., M-FOCUSS and M-SBL, were proposed [10,11]. Recently, the emerging application that exploits the sparse representation technique to achieve the DOA estimation has been increasingly attractive. Compared with the alternative methods, sparse representation-based DOA estimation methods can achieve higher angle resolution and better adapt to some challenging scenarios such as low signal-to-noise ratio (SNR), as discussed in [7]. Among the abovementioned sparse signal reconstruction methods, $l_1$-norm minimization is widely used for sparse DOA estimation because of its desirable solution accuracy and reasonable computing speed [7,10,12]. The DOA estimation algorithms using $l_1$-norm minimization include $l_1$-SRACV (sparse representation array covariance vectors) [13], real-valued (RV) $l_1$-SRACV [14], and the iterative algorithm with reweighting $l_1$-norm minimization [15]. All of these methods involve the SMV problem when recovering the signal. Aiming at the MMV case, the $l_1$-SVD (singular value decomposition) [7], RV $l_1$-SVD [16], reweighted $l_1$-SVD [17] and CMSR (covariance matrix sparse representation) [10] were proposed for DOA estimation. However, the computational complexity of the $l_1$-norm minimization is relatively large compared with that of the fast sparse representation method named smoothed $l_0$-norm algorithm, which was proposed in [18,19] under the SMV circumstance, and applied in [20,21] for the imaging and the power-line communications. To achieve fast sparse DOA estimation, by designing a reweighted continuous function in the SMV case, the reweighted smoothed $l_0$-norm (RSL0) algorithm was proposed in [22] based on the covariance vector. With enlarged array aperture and better angle estimation performance, the RSL0 algorithm [22] is about two orders of magnitude faster than the $l_1$-norm minimization based methods.

The previously mentioned DOA estimation methods are all proposed under the ideal circumstance with Gaussian white noise. However, the noise is often correlated in practical, thus the additive noise is spatially modeled as Gaussian colored noise rather than white noise. It has been verified that, with colored noise, the performance of the conventional algorithms including the methods mentioned above, is seriously degraded, except for the specially designed methods such as those based on fourth-order cumulants (FOC) [12,23,24]. These FOC-based methods are computationally expensive when they are applied to estimate the DOAs for MIMO radar. In [25], the reduced-dimension (RD) FOC-based sparse representation method (RD $l_1$-SRFOC) was proposed for DOA estimation in MIMO radar with array errors. RD $l_1$-SRFOC is an $l_1$-norm minimization based sparse representation method, which uses the high-order cumulants to deal with the colored noise and simultaneously solves the problem of mutual coupling.

It can be concluded that the existing sparse DOA estimation should improve the following problems: (I) for the $l_1$-norm minimization based methods, the computation complexities need to be lowered down and the computing time needs speeding up; (II) the fast sparse DOA estimation algorithm, i.e., the reweighted smoothed $l_0$-norm algorithm [22], is based on the covariance vector that is deemed the SMV problem. The fast sparse algorithm tailored for the MMV case needs to be studied and put forward; and (III) fast sparse algorithms for solving the colored Gaussian noise have not been available so far.

To solve the three problems stated above, a joint smoothed $l_0$-norm algorithm for DOA estimation in the MIMO radar is proposed in this paper. The new sparse algorithm first obtains a low-complexity FOC-based data matrix to eliminate the white or colored Gaussian noises. Secondly, the proposed algorithm designs a joint smoothed function tailored for the MMV case. Thirdly, with the MMV-based joint smoothed function, a joint smoothed $l_0$-norm sparse representation framework is constructed. Finally, the corresponding gradient-based sparse signal reconstruction is designed, and then DOA estimation is achieved. The proposed algorithm is a fast sparse DOA estimation algorithm, which can solve the multiple measurement vector problem, and perform well for both of the white and colored Gaussian noise environments. The proposed algorithm is about two orders of magnitude faster than the $l_1$-norm minimization based methods. This is because it approximates the $l_0$-norm by a joint smoothed function and performs the gradient-based steepest ascent scheme, which avoids the convex optimization problem involved in the $l_1$-norm minimization. The proposed algorithm can provide better DOA estimation performance than $l_1$-SVD, RV $l_1$-SVD and RV $l_1$-SRACV methods.

The rest of this paper is organized as follows. Section 2 presents the MIMO radar system and array signal models. In Section 3, the implementation process of the proposed algorithm is described in detail. Section 4 gives some related remarks and discussions regarding the advantages, the extended applications, the computational complexity and the noise elimination of the proposed algorithm. The experimental results of different methods are compared in Section 5. Finally, Section 6 concludes this paper.

Throughout the paper, the notations ${(·)}^{*}$, ${(·)}^{\mathrm{T}}$, ${(·)}^{\mathrm{H}}$, ${(·)}^{-1}$ and ${(·)}^{+}$ denote conjugation, transpose, conjugate transpose, inverse and pseudo-inverse, respectively. ${\mathbf{I}}_K$ represents a $K\times K$-dimensional unit matrix. $\left|\right|·{\left|\right|}_0$ and $\left|\right|·{\left|\right|}_2$ denote the $l_0$-norm and the $l_2$-norm, respectively. In addition, we use ⊗, $\mathrm{E}(·)$ and $\mathrm{vec}(·)$ to separately indicate the Kronecker product, the expectation and the vectorization operators.

2. Problem Formulation

A narrowband monostatic MIMO radar system is considered, and it is shown in Figure 1. The transmit and the receive arrays of this system are equipped with M transmit and N receive antennas, respectively. The arrays are both half-wavelength spaced uniform linear arrays (ULAs), whose effects of the array errors, including mutual coupling and gain-phase errors, can be ignored. The transmitting antennas transmit M orthogonal narrowband waveforms, such as BPSK (binary phase shift keying) modulated signal waveforms. In the far field, there are P uncorrelated targets regarded as point scatterers at the same range. In addition, it is assumed that $P\le M+N-2$ [22,25,26]. Thus, at the receive array, the N antennas are impinged by the echo signals reflected by the P targets. Because of the closely located arrays in the monostatic MIMO radar, for the pth target, DOA and DOD are the same and can be denoted by ${\theta }_p$. Then, by matched filtering operation, the $N\times 1$ dimensional complex envelop of the output of the mth carrier matched filter is expressed as [14,27]

$$\begin{array}{c}\hfill {\mathbf{x}}_m\left(t\right)=\sum _{p=1}^P{\mathbf{a}}_r\left({\theta }_p\right){\mathbf{a}}_{tm}^{\mathtt{T}}\left({\theta }_p\right)s_p\left(t\right)+{\mathbf{n}}_m\left(t\right),\end{array}$$

(1)

where ${\mathbf{a}}_r\left({\theta }_p\right)={\left[1,e^{j\pi sin\left({\theta }_p\right)},e^{j\pi 2sin\left({\theta }_p\right)},\dots ,e^{j\pi (N-1)sin\left({\theta }_p\right)}\right]}^{\mathtt{T}}$ is the receive steering vector, ${\mathbf{a}}_{tm}$ is the mth element of the transmit steering vector ${\mathbf{a}}_t\left({\theta }_p\right)={\left[1,e^{j\pi sin\left({\theta }_p\right)},e^{j\pi 2sin\left({\theta }_p\right)},\dots ,e^{j\pi (M-1)sin\left({\theta }_p\right)}\right]}^{\mathtt{T}}$, $s_p\left(t\right)$ contains the target reflection coefficient and the transmitted baseband signal such as non-circular signal, and ${\mathbf{n}}_m\left(t\right)$ is the noise vector after the mth matched filter. After all the matched filters, the received data vector, i.e., the vector composed of the outputs of the M matched filters, is given by [2,25,28]

$$\begin{array}{cc}\hfill \mathbf{x}\left(t\right)& =\mathrm{vec}({[{\mathbf{x}}_1^{\mathtt{T}}\left(t\right),\dots ,{\mathbf{x}}_M^{\mathtt{T}}\left(t\right)]}^{\mathtt{T}})\hfill \\ & =\mathbf{A}\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right),\hfill \end{array}$$

(2)

where $\mathbf{x}\left(t\right)\in {\mathbb{C}}^{MN\times 1}$ and $\mathbf{n}\left(t\right)={[{\mathbf{n}}_1^{\mathtt{T}}\left(t\right),\dots ,{\mathbf{n}}_M^{\mathtt{T}}\left(t\right)]}^{\mathtt{T}}\in {\mathbb{C}}^{MN\times 1}$ is the zero-mean Gaussian white or colored noise vector. $\mathbf{s}\left(t\right)={\left[s_1\left(t\right),s_2\left(t\right),\dots ,s_P\left(t\right)\right]}^{\mathtt{T}}\in {\mathbb{C}}^{P\times 1}$ is the reflected signal vector after the M carrier matched filters, and it is assumed to be statistically independent and non-Gaussian with zero-mean. The target model is considered as the classical Swerling case II, namely, the radar cross section (RCS) fluctuations are constant during a snapshot period and vary independently from snapshot to snapshot [29]. $\mathbf{n}\left(t\right)$ and $\mathbf{s}\left(t\right)$ are independent of each other. Moreover, $\mathbf{A}\in {\mathbb{C}}^{MN\times P}$ is the transmit–receive steering matrix, and its detailed expression is

$$\begin{array}{cc}\hfill \mathbf{A}& =[{\mathbf{a}}_t\left({\theta }_1\right)\otimes {\mathbf{a}}_r\left({\theta }_1\right),\dots ,{\mathbf{a}}_t\left({\theta }_P\right)\otimes {\mathbf{a}}_r\left({\theta }_P\right)]\hfill \\ & =[\mathbf{a}\left({\theta }_1\right),\dots ,\mathbf{a}\left({\theta }_P\right)],\hfill \end{array}$$

(3)

where $\mathbf{a}\left({\theta }_p\right)={\mathbf{a}}_t\left({\theta }_p\right)\otimes {\mathbf{a}}_r\left({\theta }_p\right)\in {\mathbb{C}}^{MN\times 1}$ is the transmit–receive steering vector for $p=1,2,\dots ,P$. Therefore, by collecting J snapshots, the $MN\times J$ dimensional received data matrix in monostatic MIMO radar is represented as

$$\mathbf{X}=\mathbf{A}\mathbf{S}+\mathbf{N},$$

(4)

where $\mathbf{S}=[\mathbf{s}\left(t_1\right),\dots ,\mathbf{s}\left(t_J\right)]\in {\mathbb{C}}^{P\times J}$ and $\mathbf{N}=[\mathbf{n}\left(t_1\right),\dots ,\mathbf{n}\left(t_J\right)]\in {\mathbb{C}}^{MN\times J}$ are the signal matrix and complex Gaussian noise matrix, respectively. Based on the received data $\mathbf{X}$, sparse DOA estimation can be viewed as the signal reconstruction that subjects to [14]

$$\mathbf{X}={\mathbf{A}}_{\widehat{\theta }}{\mathbf{S}}_{\widehat{\theta }}+\mathbf{N},$$

(5)

where ${\mathbf{A}}_{\widehat{\theta }}=[\mathbf{a}({\widehat{\theta }}_1),\dots ,\mathbf{a}({\widehat{\theta }}_L)]$ is the complete dictionary with the discrete sample grid ${\{{\widehat{\theta }}_i\}}_{i=1}^L$, $L\gg P$. ${\mathbf{S}}_{\widehat{\theta }}$ has the same row support with $\mathbf{S}$. Let ${\mathbf{s}}_{\widehat{\theta }}=[\mathbf{s}({\widehat{\theta }}_1),\dots ,\mathbf{s}({\widehat{\theta }}_L)]$ be a sparse vector whose ith element can be equal to the $l_2$-norm of the ith row in ${\mathbf{S}}_{\widehat{\theta }}$. To obtain the sparsest solution of Equation (5), an ideal constraint is the $l_0$-norm method by minimizing the nonzero number of ${\mathbf{s}}_{\widehat{\theta }}$, which can be expressed as

$$min\parallel {\mathbf{s}}_{\widehat{\theta }}{\parallel }_0,\mathrm{s}.\mathrm{t}.\mathbf{X}={\mathbf{A}}_{\widehat{\theta }}{\mathbf{S}}_{\widehat{\theta }}+\mathbf{N}.$$

(6)

Unfortunately, the $l_0$-norm minimization method in Equation (6) is NP-hard. To solve this problem, conventional sparse DOA estimation methods [7,10,13,14,15,17] can obtain the DOAs with the $l_1$-norm minimization. However, the computational complexity of the $l_1$-norm minimization is relatively large. In [22], by solving the array covariance vector based SMV sparse signal reconstruction problem, the reweighted smoothed $l_0$-norm algorithm greatly improves the computation speed. To develop a fast sparse DOA estimation algorithm tailored for the MMV problem existed in the MIMO radar systems, a joint smoothed $l_0$-norm algorithm is proposed in the following.

3. The Proposed Algorithm

3.1. High-Order Cumulants Based Data Matrix Construction

Based on Equation (3), the detailed expression of the steering vector $\mathbf{a}\left({\theta }_p\right)$ is

$$\begin{array}{cc}\hfill \mathbf{a}\left({\theta }_p\right)& ={[{\mathbf{a}}_r^{\mathtt{T}}\left({\theta }_p\right),e^{j\pi sin\left({\theta }_p\right)}{\mathbf{a}}_r^{\mathtt{T}}\left({\theta }_p\right),e^{j\pi (M-1)sin\left({\theta }_p\right)}{\mathbf{a}}_r^{\mathtt{T}}\left({\theta }_p\right)]}^{\mathtt{T}}\hfill \\ & ={\left[1,\dots ,z^{N-1},z,\dots ,z^N,\dots ,z^{M-1},\dots ,z^{M+N-2}\right]}^{\mathtt{T}},\hfill \end{array}$$

(7)

with $z=e^{j\pi sin\left({\theta }_p\right)}$. It can be observed that there are only $M+N-1$ distinct elements in $\mathbf{a}\left({\theta }_p\right)$; thus, the dimension of the received data vector $\mathbf{x}\left(t\right)$ can be reduced. Let $\mathbf{b}\left({\theta }_p\right)$ be a new steering vector composed of the distinct elements. It is written as

$$\mathbf{b}\left({\theta }_p\right)={[1,e^{j\pi sin\left({\theta }_p\right)},e^{j\pi 2sin\left({\theta }_p\right)},\dots ,e^{j\pi (M+N-2)sin\left({\theta }_p\right)}]}^{\mathtt{T}},$$

(8)

where $\mathbf{b}\left({\theta }_p\right)\in {\mathbb{C}}^{(M+N-1)\times 1}$. According to the element structures of $\mathbf{a}\left({\theta }_p\right)$ and $\mathbf{b}\left({\theta }_p\right)$, their relationship can be expressed as [25]

$$\mathbf{a}\left({\theta }_p\right)=\mathbf{Gb}\left({\theta }_p\right),$$

(9)

where $\mathbf{G}={[{\mathbf{L}}_0^{\mathrm{T}},{\mathbf{L}}_1^{\mathrm{T}},\cdots ,{\mathbf{L}}_{M-1}^{\mathrm{T}}]}^{\mathrm{T}}\in {\mathbb{C}}^{MN\times (M+N-1)}$ and ${\mathbf{L}}_m=[{\mathbf{0}}_{N\times m},{\mathbf{I}}_N,{\mathbf{0}}_{N\times (M-m-1)}]\in {\mathbb{C}}^{N\times (M+N-1)}$ for $m=0,1,...,M-1$. Based on the relationship in Equation (9), a reduced-dimensional matrix can be defined as [25]

$$\mathbf{R}={\left({\mathbf{G}}^{\mathrm{H}}\mathbf{G}\right)}^{(-\frac{1}{2})}{\mathbf{G}}^{\mathrm{H}}\in {\mathbb{C}}^{(M+N-1)\times MN}.$$

(10)

Therefore, the dimensional reduction for the received data is carried out as follows [25]:

$$\begin{array}{cc}\hfill \tilde{\mathbf{x}}\left(t\right)& =\mathbf{R}\mathbf{G}\mathbf{B}\mathbf{s}\left(t\right)+\mathbf{R}\mathbf{n}\left(t\right)\hfill \\ & ={\left({\mathbf{G}}^{\mathrm{H}}\mathbf{G}\right)}^{(-\frac{1}{2})}{\mathbf{G}}^{\mathrm{H}}\mathbf{G}\mathbf{B}\mathbf{s}\left(t\right)+\tilde{\mathbf{n}}\left(t\right)\hfill \\ & =\mathbf{F}\mathbf{B}\mathbf{s}\left(t\right)+\tilde{\mathbf{n}}\left(t\right),\hfill \end{array}$$

(11)

where $\tilde{\mathbf{x}}\left(t\right)=\mathbf{R}\mathbf{x}\left(t\right)\in {\mathbb{C}}^{(M+N-1)\times 1}$, $\mathbf{B}=[\mathbf{b}\left({\theta }_1\right),\dots ,\mathbf{b}\left({\theta }_P\right)]\in {\mathbb{C}}^{(M+N-1)\times P}$ and $\tilde{\mathbf{n}}\left(t\right)\in {\mathbb{C}}^{(M+N-1)\times 1}$ are the new low-dimensional received data vector, transmit–receive steering matrix and noise vector, respectively. In addition, $\mathbf{F}={\left({\mathbf{G}}^{\mathrm{H}}\mathbf{G}\right)}^{\frac{1}{2}}$ is a known diagonal matrix, and its $(i,i)$th element can be expressed as

$$\begin{array}{c}\hfill \mathbf{F}(i,i)=\left\{\begin{array}{c}\sqrt{i}i=1,2,\dots ,\beta \hfill \\ \sqrt{\beta }i=\beta +1,\beta +2,\dots ,\beta +\tau -1\hfill \\ \sqrt{M+N-i}i=\beta +\tau ,\dots ,M+N-1,\hfill \end{array}\right.\end{array}$$

(12)

with $\beta =\min (M,N)$ and $\tau =|M-N|+2$. Note that the noise vector $\tilde{\mathbf{n}}\left(t\right)$ remains Gaussian after the dimensional reduction, for asymptotic normal distribution has the invariance speciality of linear transformation. To eliminate the Gaussian white noise or colored noise in $\tilde{\mathbf{x}}\left(t\right)$, fourth-order cumulant is exploited under the circumstance of collecting J snapshots. The definition of the FOC is given as follows [12,24,25]

$$\begin{array}{cc}\hfill C_{4\tilde{x}}(k_1,k_2,k_3,k_4)& =\mathrm{cum}\left\{{\tilde{\mathbf{x}}}_{k_1},{\tilde{\mathbf{x}}}_{k_2}^{*},{\tilde{\mathbf{x}}}_{k_3},{\tilde{\mathbf{x}}}_{k_4}^{*}\right\}=\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_1}{\tilde{\mathbf{x}}}_{k_2}^{*}{\tilde{\mathbf{x}}}_{k_3}{\tilde{\mathbf{x}}}_{k_4}^{*}\right)-\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_1}{\tilde{\mathbf{x}}}_{k_2}^{*}\right)\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_3}{\tilde{\mathbf{x}}}_{k_4}^{*}\right)\hfill \\ & -\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_1}{\tilde{\mathbf{x}}}_{k_3}\right)\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_2}^{*}{\tilde{\mathbf{x}}}_{k_4}^{*}\right)-\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_1}{\tilde{\mathbf{x}}}_{k_4}^{*}\right)\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_2}^{*}{\tilde{\mathbf{x}}}_{k_3}\right),\hfill \end{array}$$

(13)

where $\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_1}{\tilde{\mathbf{x}}}_{k_2}^{*}{\tilde{\mathbf{x}}}_{k_3}{\tilde{\mathbf{x}}}_{k_4}^{*}\right)\approx \frac{1}{J}{\sum }_{t=1}^J{\tilde{\mathbf{x}}}_{k_1}\left(t\right){\tilde{\mathbf{x}}}_{k_2}^{*}\left(t\right){\tilde{\mathbf{x}}}_{k_3}\left(t\right){\tilde{\mathbf{x}}}_{k_4}^{*}\left(t\right)$, $\mathrm{E}\left({\tilde{\mathbf{x}}}_{k_1}{\tilde{\mathbf{x}}}_{k_2}^{*}\right)\approx \frac{1}{J}{\sum }_{t=1}^J{\tilde{\mathbf{x}}}_{k_1}\left(t\right){\tilde{\mathbf{x}}}_{k_2}^{*}\left(t\right)$, ${\tilde{\mathbf{x}}}_{k_i}$ stands for the $k_i$th element in $\tilde{\mathbf{x}}$, and $1\le k_i\le M+N-1$ for $i=1,2,3,4$. With the assumptions about the noise in Equation (2), the fourth-order cumulant of Gaussian noise has the following property:

$$C_{4\tilde{n}}(k_1,k_2,k_3,k_4)=\mathrm{cum}\left\{{\tilde{\mathbf{n}}}_{k_1},{\tilde{\mathbf{n}}}_{k_2}^{*},{\tilde{\mathbf{n}}}_{k_3},{\tilde{\mathbf{n}}}_{k_4}^{*}\right\}=0,$$

(14)

where $C_{4\tilde{n}}(k_1,k_2,k_3,k_4)$ represents the fourth-order cumulant of $\tilde{\mathbf{n}}$ corresponding to the indices $(k_1,k_2,k_3,k_4)$. For the signal assumed in Equation (2), the fourth-order cumulant satisfies

$$\begin{array}{c}\hfill C_{4s}(p_1,p_2,p_3,p_4)=\left\{\begin{array}{c}{c}_p{p}_1=p_2=p_3=p_4=p\hfill \\ 0otherwise,\hfill \end{array}\right.\end{array}$$

(15)

where $c_p=\mathrm{cum}\{{\mathbf{s}}_p,{\mathbf{s}}_p^{*},{\mathbf{s}}_p,{\mathbf{s}}_p^{*}\}$ is the FOC of the signal for the pth target, and $1\le p_i\le P$ for $i=1,2,3,4$. Based on Equations (14) and (15), an FOC-based data matrix can be constructed with its $(k_1,k_2)$th element being obtained from

$$\begin{array}{cc}\hfill \mathbf{Y}(k_1,k_2)& =C_{4\tilde{x}}(k_1,k_2,k_2,k_2)\hfill \\ & =\mathrm{cum}\left\{{\tilde{\mathbf{x}}}_{k_1},{\tilde{\mathbf{x}}}_{k_2}^{*},{\tilde{\mathbf{x}}}_{k_2},{\tilde{\mathbf{x}}}_{k_2}^{*}\right\}\hfill \\ & =\sum _{p=1}^P{\left[\mathbf{F}\mathbf{b}\left({\theta }_p\right)\right]}_{k_1}{\left[\mathbf{F}{\mathbf{b}}^{*}\left({\theta }_p\right)\right]}_{k_2}{\left|\mathbf{F}(k_2,k_2)\right|}^2{c}_p,\hfill \end{array}$$

(16)

where ${\left[\mathbf{F}\mathbf{b}\left({\theta }_p\right)\right]}_{k_1}$ and ${\left[\mathbf{F}{\mathbf{b}}^{*}\left({\theta }_p\right)\right]}_{k_2}$ are the $k_1$th and the $k_2$th elements in $\mathbf{F}\mathbf{b}\left({\theta }_p\right)$ and $\mathbf{F}{\mathbf{b}}^{*}\left({\theta }_p\right)$. As a result, the data matrix $\mathbf{Y}$ can be expressed as

$$\mathbf{Y}=\mathbf{F}\mathbf{B}{\mathbf{D}}_s{\left({\mathbf{F}}^3\mathbf{B}\right)}^{\mathtt{H}},$$

(17)

where ${\mathbf{D}}_s=\mathrm{diag}(c_1,c_2,\dots ,c_P)$ is a diagonal matrix composed of the fourth-order cumulants for ${\mathbf{s}}_p$, $p=1,2,\dots ,P$. Then, the singular value decomposition of $\mathbf{Y}$ can be performed as follows:

$$\mathbf{Y}=\mathbf{U}\mathsf{\Lambda }{\mathbf{V}}^{\mathtt{H}},$$

(18)

where $\mathbf{U}$ is the left singular vector, $\mathbf{V}$ is the right singular vector, and $\mathsf{\Lambda }=\mathrm{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{M+N-1})$ is the singular value matrix with ${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_{M+N-1}$. Therefore, the signal subspace ${\mathbf{V}}_s$ is obtained by extracting the first P right singular vectors that correspond to $({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_P)$. In addition, the last $M+N-1-P$ left singular vectors corresponding to $({\lambda }_{P+1},{\lambda }_{P+2},\dots ,{\lambda }_{M+N-1})$ are used to make up the noise subspace ${\mathbf{U}}_n\in {\mathbb{C}}^{(M+N-1)\times (M+N-1-P)}$. With ${\mathbf{V}}_s\in {\mathbb{C}}^{(M+N-1)\times P}$, $\mathbf{Y}$ can be transformed into a dimension-reduced data matrix, that is

$$\begin{array}{c}\hfill {\mathbf{Y}}_s=\mathbf{Y}{\mathbf{V}}_s=\mathbf{F}\mathbf{B}{\mathbf{D}}_s{\left({\mathbf{F}}^3\mathbf{B}\right)}^{\mathtt{H}}{\mathbf{V}}_s,\end{array}$$

(19)

where ${\mathbf{Y}}_s\in {\mathbb{C}}^{(M+N-1)\times P}$.

3.2. Designs of Joint Smoothed Function and Joint Smoothed $l_0$-Norm Framework for MMV Case

In the sparse representation framework for DOA estimation, a complete dictionary containing all interest DOAs is required. Thus, let ${\{{\widehat{\theta }}_i\}}_{i=1}^L$ be the discrete sample grid. According to Equation (19), the complete dictionary can be constructed as ${\mathbf{B}}_{\widehat{\theta }}=[\mathbf{b}({\widehat{\theta }}_1),\mathbf{b}({\widehat{\theta }}_2),\dots ,\mathbf{b}({\widehat{\theta }}_L)]\in {\mathbb{C}}^{(M+N-1)\times L}$. The DOA estimation is considered as a sparse signal reconstruction problem that is subject to

$${\mathbf{Y}}_s=\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}{\mathbf{T}}_{\widehat{\theta }},$$

(20)

where ${\mathbf{T}}_{\widehat{\theta }}\in {\mathbb{C}}^{L\times P}$ is a sparse matrix and has the same row support with $\mathbf{T}={\mathbf{D}}_s{\left({\mathbf{F}}^3\mathbf{B}\right)}^{\mathtt{H}}{\mathbf{V}}_s$. Then, DOA estimation can be achieved by measuring the sparsity of ${\mathbf{T}}_{\widehat{\theta }}$. Note that, in the sparse signal recovery problem in Equation (20), the low-dimensional data ${\mathbf{Y}}_s$ is a matrix rather than a vector. Therefore, recovering ${\mathbf{T}}_{\widehat{\theta }}$ involves the MMV problem. To develop a fast sparse DOA estimation algorithm that can solve the MMV problem, in the following, the proposed joint smoothed $l_0$-norm algorithm designs a new joint smoothed function and then derives its new steepest ascent scheme. For obtaining the smoothed estimation of the sparsity of ${\mathbf{T}}_{\widehat{\theta }}$, we first design a joint continuous function as follows

$$f_{\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]=\exp \{-[\mid {\mathbf{T}}_{\widehat{\theta }}(i,1){\mid }^2+\mid {\mathbf{T}}_{\widehat{\theta }}(i,2){\mid }^2+\dots +\mid {\mathbf{T}}_{\widehat{\theta }}(i,P){\mid }^2]/2P{\sigma }^2\},$$

(21)

where ${\mathbf{T}}_{\widehat{\theta }}(i,p)$ is the $(i,p)$th element of ${\mathbf{T}}_{\widehat{\theta }}$ for $p=1,2,\dots ,P$ and $i=1,2,\dots ,L$. Thus, the result of the designed joint continuous function $f_{\sigma }[{\mathbf{T}}_{\widehat{\theta }}(i,:)]$ can be expressed as follows:

$$\underset{\sigma \to 0}{\lim }{f}_{\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]=\left\{\begin{array}{c}1{\mathbf{T}}_{\widehat{\theta }}^{l_2}\left(i\right)=0\\ 0,{\mathbf{T}}_{\widehat{\theta }}^{l_2}\left(i\right)\ne 0,\end{array}\right.$$

(22)

where ${\mathbf{T}}_{\widehat{\theta }}^{\left(l_2\right)}\in {\mathbb{C}}^{L\times 1}$ is the $L\times 1$ dimensional vector, and $[{\mathbf{T}}_{\widehat{\theta }}^{\left(l_2\right)}]\left(i\right)=\parallel {\mathbf{T}}_{\widehat{\theta }}(i,:){\parallel }_2$. Namely, the ith element of ${\mathbf{T}}_{\widehat{\theta }}^{\left(l_2\right)}$ is the $l_2$-norm of the ith row in ${\mathbf{T}}_{\widehat{\theta }}$. Then, a series of weights $[r_{w1},r_{w2},\dots ,r_{wL}]$ are introduced into the joint continuous function $f_{\sigma }[{\mathbf{T}}_{\widehat{\theta }}(i,:)]$. Therefore, in the proposed joint smoothed $l_0$-norm algorithm, the final design of the joint smoothed function is

$$\begin{array}{cc}\hfill f_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]& =\exp \{-\frac{r_{wi}}{2P{\sigma }^2}\sum _{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]\}\hfill \\ & =e^{-\frac{r_{wi}}{2P{\sigma }^2}{\sum }_{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]},\hfill \end{array}$$

(23)

where $r_{wi}$ is the weight coefficient, which is obtained from $r_{wi}=\parallel {[\mathbf{F}\mathbf{b}({\widehat{\theta }}_i)]}^{\mathtt{H}}{\mathbf{U}}_n{\parallel }_2/r_{max}$ with $r_{max}=\max \{\parallel {[\mathbf{F}\mathbf{b}({\widehat{\theta }}_1)]}^{\mathtt{H}}{\mathbf{U}}_n{\parallel }_2,\parallel {[\mathbf{F}\mathbf{b}({\widehat{\theta }}_2)]}^{\mathtt{H}}{\mathbf{U}}_n{\parallel }_2,\dots ,\parallel {[\mathbf{F}\mathbf{b}({\widehat{\theta }}_L)]}^{\mathtt{H}}{\mathbf{U}}_n{\parallel }_2\}$. In Equation (23), the preset parameter $\sigma$ is known and adjusts the smoothness of $f_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]$. In addition, P stands for the total number of the columns in the data matrix ${\mathbf{Y}}_s$, namely, the column number of the sparse matrix ${\mathbf{T}}_{\widehat{\theta }}$. For the true target DOA ${\theta }_p$, $0<r_{wp}\ll 1$ because of the orthogonality between the signal subspace and the noise subspace [5]. Therefore, the approximate value of the joint smoothed function $f_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]$ corresponding to ${\widehat{\theta }}_i$ can be calculated as follows:

$$f_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]\approx \left\{\begin{array}{c}1r_{wi}\to 1,{\mathbf{T}}_{\widehat{\theta }}^{l_2}\left(i\right)\ll \sigma \\ 0r_{wi}\ll 1,{\mathbf{T}}_{\widehat{\theta }}^{l_2}\left(i\right)\gg \sigma ,\end{array}\right.$$

(24)

where a small $\sigma$ is required to guarantee the approximation. Based on Equation (24), the sparsity of ${\mathbf{T}}_{\widehat{\theta }}$ can be measured by

$$\begin{array}{c}\hfill \left|\right|{\mathbf{T}}_{\widehat{\theta }}^{l_2}{\left|\right|}_0\approx L-\sum _{i=1}^L{f}_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]=L-F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2}),\end{array}$$

(25)

where $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})={\sum }_{i=1}^L{f}_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]$. According to Equation (25), to recover the sparse matrix ${\mathbf{T}}_{\widehat{\theta }}$, the minimization of $\left|\right|{\mathbf{T}}_{\widehat{\theta }}^{l_2}{\left|\right|}_0$ can be turned into the maximization of $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$. As a consequence, the MMV-based joint smoothed $l_0$-norm sparse representation framework is represented as

$$\underset{{\mathbf{T}}_{\widehat{\theta }}}{max}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2}),\mathrm{s}.\mathrm{t}.{\mathbf{Y}}_s=\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}{\mathbf{T}}_{\widehat{\theta }}.$$

(26)

3.3. Design of the MMV-Based Signal Reconstruction in the Joint Smoothed $l_0$-Norm Algorithm

For the joint smoothed function $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$, the steepest ascent based on the gradient is designed to recover the sparse signal. Aiming at the sparse signal reconstruction in Equation (26) that involves the MMV problem, an initial matrix is needed for the joint smoothed $l_0$-norm algorithm to start the steepest ascent process. The initial matrix ${\mathbf{U}}_0^{\widehat{\theta }}$ can be chosen as the solution of Equation (26) with $\sigma \to \infty$, i.e., the minimum $l_2$-norm solution calculated by

$${\mathbf{U}}_0^{\widehat{\theta }}={\left(\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}\right)}^{+}{\mathbf{Y}}_s\in {\mathbb{C}}^{(M+N-1)\times P},$$

(27)

where ${\left(\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}\right)}^{+}$ is the pseudo-inverse of $\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}$. When recovering the sparse signal for Equation (26), the selected value of $\sigma$ has an effect on the computation time and the solution accuracy. The larger $\sigma$ is, the smoother $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ will be. Larger $\sigma$ can result in less local maxima and makes it more easy to maximize $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$. In this case, the computation process can be speeded up. However, large $\sigma$ causes the deviation of the approximation in Equation (25), and it will further decrease the precision of the signal reconstruction. To guarantee both of the speed and the accuracy, a decreasing sequence for $\sigma$, i.e., $[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_K]$ with ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_K$, which is known as the graduated non-convexity (GNC) approach [18,22], can be exploited. Even though large $\sigma$ is desired to avoid a lot of local maxima and improve the speed, the value of ${\sigma }_1$ in the decreasing $\sigma$ sequence may be set as one to four times of the maximum absolute value of ${\mathbf{U}}_0^{\widehat{\theta }}$. This is because in $f_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]$, when $\sigma >4\max \{|{\mathbf{U}}_0^{\widehat{\theta }}|\}$, $\sigma$ acts virtually like infinity for all elements with $e^{-\frac{1}{2P{\sigma }^2}{\sum }_{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]}>0.96\approx 1$. Then, the other values in the $\sigma$ sequence can be chosen as ${\sigma }_{k+1}=\alpha {\sigma }_k$ for $k=1,2,\dots ,K-1$. To make the steepest ascent algorithm be not far from the actual maximum, and be less likely to get trapped into local maxima, $\sigma$ is required to change slowly. This guarantees that the GNC approach can work. In general, the decreasing factor $\alpha$ satisfies $0.5\le \alpha <1$. ${\sigma }_K$ should be a small value determined by the desired accuracy, and a parameter ${\sigma }_{off}$ can be used to set the lower limit of the decreasing $\sigma$ sequence.

For each ${\sigma }_k$ in the decreasing $\sigma$ sequence, graduated non-convexity based Q iterations are carried out to maximize $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$. In each iteration, ${\mathbf{T}}_{\widehat{\theta }}$ is first updated as the following form:

$${\mathbf{T}}_{\widehat{\theta }}\leftarrow {\mathbf{T}}_{\widehat{\theta }}+\tilde{\mu }\nabla F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})={\mathbf{T}}_{\widehat{\theta }}-\mu \Delta {\mathbf{T}}_{\widehat{\theta }},$$

(28)

where $\tilde{\mu }=\mu {\sigma }^2$ is a step-size parameter, and $\mu$ should be $\mu \ge 1$. In addition, $\nabla F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ is the gradient of $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ corresponding to ${\mathbf{T}}_{\widehat{\theta }}$. In the following, $\Delta {\mathbf{T}}_{\widehat{\theta }}$ will be derived in detail. Let the detailed expression of ${\mathbf{T}}_{\widehat{\theta }}$ be

$${\mathbf{T}}_{\widehat{\theta }}={[{\mathbf{t}}_1^{\mathtt{T}},{\mathbf{t}}_2^{\mathtt{T}},\dots ,{\mathbf{t}}_L^{\mathtt{T}}]}^{\mathtt{T}},$$

(29)

where ${\mathbf{t}}_i\in {\mathbb{C}}^{1\times P}$ is the ith row of ${\mathbf{T}}_{\widehat{\theta }}$. According to Equation (28), ${\mathbf{t}}_i$ can be written as

$${\mathbf{t}}_i\leftarrow {\mathbf{t}}_i+\tilde{\mu }{\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})={\mathbf{t}}_i-\mu \Delta {\mathbf{t}}_i,$$

(30)

where ${\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ stands for the $1\times P$ dimensional gradient vector of $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ in regard to ${\mathbf{t}}_i=[{\mathbf{T}}_{\widehat{\theta }}(i,1),{\mathbf{T}}_{\widehat{\theta }}(i,2),\dots ,{\mathbf{T}}_{\widehat{\theta }}(i,P)]$. The pth element of ${\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ is detailedly expressed as

$$\begin{array}{cc}\hfill {\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})\left(p\right)& =\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial {\mathbf{t}}_{ip}\hfill \\ & =\partial \sum _{i=1}^L{e}^{-\frac{r_{wi}}{2P{\sigma }^2}{\sum }_{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]}/\partial {\mathbf{T}}_{\widehat{\theta }}(i,p),\hfill \end{array}$$

(31)

where ${\mathbf{t}}_{ip}=\tilde{a}+\tilde{b}j$ is the pth parameter of ${\mathbf{t}}_i$ with the real component $\tilde{a}$ and imaginary component $\tilde{b}$. $\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial {\mathbf{t}}_{ip}=[\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial \tilde{a}]+[\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial \tilde{b}]j$ denotes the 1st-order derivative of $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$. As a result, ${\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})\left(p\right)$ is derived as follows:

$$\begin{array}{cc}\hfill {\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})\left(p\right)& =e^{-\frac{r_{wi}}{2P{\sigma }^2}{\sum }_{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]}\partial \{-\frac{r_{wi}}{2P{\sigma }^2}\sum _{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]\}/\partial {\mathbf{T}}_{\widehat{\theta }}(i,p)\hfill \\ & =-\frac{r_{wi}{\mathbf{T}}_{\widehat{\theta }}(i,p)}{P{\sigma }^2}\exp \{-\frac{r_{wi}}{2P{\sigma }^2}\sum _{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]\}.\hfill \end{array}$$

(32)

Combining Equation (30) with Equation (32), the pth element of $\Delta {\mathbf{t}}_i$ can be obtained, that is,

$$\begin{array}{cc}\hfill \Delta {\mathbf{t}}_{ip}& =-{\sigma }^2{\nabla }_{{\mathbf{t}}_i}{F}_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})\left(p\right)\hfill \\ & =\frac{r_{wi}{\mathbf{T}}_{\widehat{\theta }}(i,p)}{P}\exp \{-\frac{r_{wi}}{2P{\sigma }^2}\sum _{p=1}^P[\mid {\mathbf{T}}_{\widehat{\theta }}(i,p){\mid }^2]\}.\hfill \end{array}$$

(33)

With $\Delta {\mathbf{t}}_{ip}$, the $1\times P$ dimensional vector $\Delta {\mathbf{t}}_i$ is written as

$$\begin{array}{cc}\hfill \Delta {\mathbf{t}}_i& =-{\sigma }^2[\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial {\mathbf{T}}_{\widehat{\theta }}(i,1),\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial {\mathbf{T}}_{\widehat{\theta }}(i,2),\dots ,\partial F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})/\partial {\mathbf{T}}_{\widehat{\theta }}(i,P)]\hfill \\ & =[\Delta {\mathbf{t}}_{i1},\Delta {\mathbf{t}}_{i2},\dots ,\Delta {\mathbf{t}}_{iP}]\hfill \end{array}$$

(34)

for $i=1,2,\dots ,L$. As a result, based on Equations (33) and (34), the sparse matrix to be recovered in each iteration is first updated as

$${\mathbf{T}}_{\widehat{\theta }}\leftarrow {\mathbf{T}}_{\widehat{\theta }}-\mu \Delta {\mathbf{T}}_{\widehat{\theta }}\in {\mathbb{C}}^{L\times P},$$

(35)

with

$$\Delta {\mathbf{T}}_{\widehat{\theta }}={[\Delta {\mathbf{t}}_1^{\mathtt{T}},\Delta {\mathbf{t}}_2^{\mathtt{T}},\dots ,\Delta {\mathbf{t}}_P^{\mathtt{T}}]}^{\mathtt{T}},$$

(36)

where $\Delta {\mathbf{T}}_{\widehat{\theta }}\in {\mathbb{C}}^{L\times P}$. Then, this iteration projects ${\mathbf{T}}_{\widehat{\theta }}$ back to the feasible set ${\mathcal{T}}_{\widehat{\theta }}=\{{\mathbf{T}}_{\widehat{\theta }}\mid {\mathbf{Y}}_s=\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}{\mathbf{T}}_{\widehat{\theta }}\},$ which satisfies Equation (20), that is,

$${\mathbf{T}}_{\widehat{\theta }}\leftarrow {\mathbf{T}}_{\widehat{\theta }}-{\left(\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}\right)}^{+}[\left(\mathbf{F}{\mathbf{B}}_{\widehat{\theta }}\right){\mathbf{T}}_{\widehat{\theta }}-{\mathbf{Y}}_s].$$

(37)

Based on Equations (27)–(37), after Q iterations for $\sigma ={\sigma }_k$, the MMV based sparse signal reconstruction in Equation (26) starts the maximization process of the joint smoothed function $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ for $\sigma ={\sigma }_{k+1}$. Let ${\mathbf{T}}_f^{\widehat{\theta }}$ be the output of ${\mathbf{T}}_{\widehat{\theta }}$ by maximizing $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$ for $\sigma ={\sigma }_K$. Consequently, the final sparse matrix is recovered as ${\mathbf{T}}_f^{\widehat{\theta }}$. In order to evaluate the reconstruction results of each row in ${\mathbf{T}}_f^{\widehat{\theta }}$, and then achieve the DOA estimation, let

$${\overline{\ell }}_{\mathbf{p}}={({\mathbf{T}}_f^{\widehat{\theta }})}^{l_2},$$

(38)

with $[{({\mathbf{T}}_f^{\widehat{\theta }})}^{l_2}]\left(i\right)=\parallel {\mathbf{T}}_f^{\widehat{\theta }}(i,:){\parallel }_2$. Therefore, by searching the spectrum of ${\overline{\ell }}_{\mathbf{p}}$, the target DOAs are obtained. The summary of the proposed joint smoothed $l_0$-norm algorithm is listed below:

Step 1:

Compute the fourth-order cumulants based matrix ${\mathbf{Y}}_s$ from Equations (11), (16) and (19).

Step 2:

Design the joint smoothed function $f_{w\sigma }\left[{\mathbf{T}}_{\widehat{\theta }}(i,:)\right]$ tailored for the MMV case as Equation (23).

Step 3:

Construct the joint smoothed $l_0$-norm sparse representation framework in Equation (26).

Step 4:

Execute the fast MMV-based sparse signal reconstruction with Equations (27)–(37).

Step 5:

Attain the DOA estimation based on Equation (38).

4. Related Remarks

Remark 1.

With the joint smoothed function $F_{w\sigma }({\mathit{T}}_{\widehat{\theta }}^{l_2})$, the proposed joint smoothed $l_0$-norm minimization based sparse signal reconstruction algorithm is a new fast sparse DOA estimation algorithm, which is tailored for the MMV problem. By designing a specific joint continuous function that exploits all data information in the data matrix ${\mathit{Y}}_s$, and then, deriving its steepest ascent process to achieve the maximization, the sparse solutions can be obtained. In the extended and improved smoothed $l_0$-norm algorithms, the proposed joint continuous function $f_{\sigma }\left[{\mathit{T}}_{\widehat{\theta }}(i,:)\right]$ in Equation (21) is applicable for all MMV cases with the form ${\mathit{Y}}_s={\widehat{\mathbf{B}}}_{\widehat{\theta }}{\mathit{T}}_{\widehat{\theta }}$, where ${\widehat{\mathbf{B}}}_{\widehat{\theta }}$ is the complete dictionary and ${\mathit{T}}_{\widehat{\theta }}$ is the $L\times P$ dimensional sparse matrix that needs to be reconstructed. Under this circumstance, gradient-based maximizing ${\sum }_{i=1}^L{f}_{\sigma }\left[{\mathit{T}}_{\widehat{\theta }}(i,:)\right]$ is similar to the derivation of Equations (27)–(37), and its $\Delta {\mathit{t}}_{ip}$ corresponding to ${\sum }_{i=1}^L{f}_{\sigma }\left[{\mathit{T}}_{\widehat{\theta }}(i,:)\right]$ is concluded as

$$\begin{array}{c}\hfill \Delta {\mathit{t}}_{ip}=\frac{{\mathit{T}}_{\widehat{\theta }}(i,p)}{P}\exp \{-\frac{1}{2P{\sigma }^2}\sum _{p=1}^P[\mid {\mathit{T}}_{\widehat{\theta }}(i,p){\mid }^2]\}\end{array}$$

(39)

for $i=1,2,\dots ,L$. With $\Delta {\mathit{t}}_{ip}$ in Equation (39), the proposed joint smoothed $l_0$-norm algorithm regarding the joint continuous function $f_{\sigma }\left[{\mathit{T}}_{\widehat{\theta }}(i,:)\right]$ can be applied to other fast sparse signal reconstruction problems that involve multiple measurement vectors.

Remark 2.

The proposed method has low computational complexity because its sparse signal recovery is a graduated non-convexity procedure, rather than the convex optimization problem involved in the $l_1$-norm minimization. When recovering the sparse signal matrix ${\mathit{T}}_{\widehat{\theta }}$, the main computational burden is caused by the $QK$ iterations, which require about $O[PLQK+(M+N-1\left)PLQK\right]$, where L is the number of the sample grid cells in the complete dictionary ${\mathbf{B}}_{\widehat{\theta }}$. For the SMV based reweighted smoothed $l_0$-norm algorithm in [22], recovering the sparse signal vector requires about $O[{(M+N-1)}^{2}LQK]$. However, $l_1$-norm minimization based DOA estimation methods, such as $l_1$-SVD and RV $l_1$-SVD, require $O[{\left(PL\right)}^3]$ and $O[\frac{1}{4}{\left(PL\right)}^3]$, respectively. In addition, the RV $l_1$-SRACV needs $O\left[\frac{1}{4}{L}^3\right]$. As L is much larger than P, $M+N-1$, Q and K, the computation speed of the proposed joint smoothed $l_0$-norm algorithm is much faster than the $l_1$-norm minimization based methods.

Remark 3.

The proposed method can deal with both of the white and colored Gaussian noises, thereby making it desirable for use in practical applications. This is because, in the proposed joint smoothed $l_0$-norm algorithm, high-order cumulants are utilized to process the reduced-dimensional data $\mathit{x}\left(t\right)$. Note that, for the dimension reduction in Equation (11), the reduced-dimensional matrix $\mathit{R}$ satisfies $\mathit{R}{\mathit{R}}^{{H}}={\mathit{I}}_{M+N-1}$. Hence, if the noise vector $\tilde{\mathit{n}}\left(t\right)$ is white Gaussian with covariance matrix ${\rho }^2{\mathit{I}}_{MN}$, it remains white Gaussian with covariance matrix $\mathit{R}{\rho }^2{\mathit{I}}_{MN}{\mathit{R}}^{{H}}={\rho }^2{\mathit{I}}_{M+N-1}$ after the transformation [2,25]. On the other hand, if the noise vector $\mathit{n}\left(t\right)$ is colored Gaussian with covariance matrix ${\mathit{V}}_a$, its covariance matrix turns into $\mathit{R}{\mathit{V}}_a{\mathit{R}}^{{H}}$ after the dimensional reduction. However, in this case, $\tilde{\mathit{n}}\left(t\right)$ is still Gaussian according to the invariance speciality of linear transformation in asymptotic normal distribution. Therefore, even though the dimensional reduction makes the change of the noise happen, the high-order cumulant process in Equation (16) can eliminate it.

5. Simulation Results

In this part, some simulations are implemented to evaluate the performance and the computation speed of the proposed method. The simulation results of the proposed method are compared with those of the $l_1$-SVD [7], RV $l_1$-SVD [16], RV $l_1$-SRACV [14] and reweighted smoothed $l_0$-norm (RSL0) [22] algorithms. In all of the simulations, a narrowband monostatic MIMO radar system with ULAs is considered. The antennas of the transmit array and the receive array are both half-wavelength spaced, and their numbers are $M=6$ and $N=6$, respectively. The transmitted signals are assumed to be BPSK modulated [27,30], moreover, they are mutually orthogonal and different from each other. The noise is white Gaussian with the covariance matrix $\mathbf{R}={\mathbf{I}}_{MN}$, or colored Gaussian with the elements in the covariance matrix being $\mathbf{R}(k_1,k_2)=0.{75}^{|k_1-k_2|}{e}^{j\pi (k_1-k_2)/2}$. By means of the minimum description length (MDL) principle or the Akaike information criterion (AIC) principle [31], the prior knowledge of the target number is assumed to be known as P.

When constructing the complete dictionary, the discrete sample grid is set from $-{90}^{\circ }$ to ${90}^{\circ }$ with the uniform interval distance $\Delta \widehat{\theta }=0.{05}^{\circ }$ for the proposed algorithm, as well as the $l_1$-SVD, RV $l_1$-SVD, RV $l_1$-SRACV and RSL0 algorithms. The parameters ${\sigma }_1$, $\mu$ and ${\sigma }_{off}$ are fixedly set at ${\sigma }_1=4\max \{|{\mathbf{U}}_0^{\widehat{\theta }}|\}$, $\mu =2.5$ and ${\sigma }_{off}=0.0007$ for the proposed algorithm and the RSL0 algorithm. The effects of the other parameters, i.e., the decreasing factor $\alpha$ and the iteration number Q, will be evaluated and given in detail in the following experiments. The definitions for the signal-to-noise ratio (SNR) and the root mean square error (RMSE) are shown, respectively, as follows:

$$\begin{array}{cc}\hfill \mathrm{SNR}& =10{\log }_{10}{\left(\right|\left|\mathbf{A}\mathbf{S}\right||}_F^2/{\left|\right|\mathbf{N}\left|\right|}_F^2),\hfill \\ \hfill \mathrm{RMSE}& =\sqrt{\frac{1}{500P}{\sum }_{i=1}^{500}{\sum }_{p=1}^P{({\widehat{\theta }}_{p,i}-{\theta }_p)}^2},\hfill \end{array}$$

(40)

where ${\widehat{\theta }}_{p,i}$ is the angle estimation of the pth target DOA ${\theta }_p$ for the ith Monte Carlo trial.

Figure 2 depicts the spatial spectrum of the proposed method for white noise and colored noise, where SNR = 5 dB and SNR = −10 dB are considered. There are three targets located at ${\theta }_1=-8.2^{\circ }$, ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$. The spatial spectrum is plotted by computing $10{\log }_{10}\left[\right|{\overline{\ell }}_{\mathbf{p}}|/\max (|{\overline{\ell }}_{\mathbf{p}}\left|\right)]$ with ${\overline{\ell }}_{\mathbf{p}}$ being the final sparse solution. Figure 2 shows that the proposed joint smoothed $l_0$-norm minimization algorithm is effective for the MMV case of DOA estimation in monostatic MIMO radar. Moreover, it is suitable for the practical radar systems with colored Gaussian noise.

Figure 3a,b evaluates the effects of the parameters $\alpha$ and Q on the DOA estimation performance of the proposed method, with $Q=3$ and $\alpha =0.5$, respectively. Moreover, $J=1000$. Three uncorrelated targets with the DOAs ${\theta }_1=-{8.2}^{\circ }$, ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$ are considered. In Figure 3, it can be verified that the proposed method keeps the same angle estimation performance when the parameters $\alpha$ and Q vary from $\alpha =0.5$ to $\alpha =0.9$ and $Q=3$ to $Q=35$. In Figure 3b, the RMSE decreases with the increase of Q when $Q\le 3$, and it reaches an asymptotic value at $Q=3$. For the iterations of the steepest ascent algorithm, we just need a region near the maximum value of $F_{w\sigma }({\mathbf{T}}_{\widehat{\theta }}^{l_2})$; thus, a fixed small positive integer $Q=3$ is applicable.

With respect to Figure 3, the effects of the parameters $\alpha$ and Q on the average computation time of the proposed method are shown in

Table 1. Average computation time for different values of the parameter $\alpha$.

$\mathit{\alpha }$	0.5	0.6	0.7	0.8	0.9
Average Computation Time (s)	0.0401	0.0498	0.0752	0.1137	0.2511

and

Table 2. Average computation time for different values of the parameter Q.

Q	3	11	19	27	35
Average Computation Time (s)	0.0399	0.1382	0.2576	0.3696	0.4738

. The simulation conditions keep the same with those of Figure 3. As the implementing processes of the analyzed methods are not affected by the type of the noise, the computation time is the same for the white and the colored Gaussian noises. It can be observed that the computation speed becomes slower when $\alpha$ and Q increase. However, larger $\alpha$ and Q are not able to improve the estimation performance of the proposed method, as shown in Figure 3a,b. As a result, small values of $\alpha$ and Q are enough to guarantee the estimation accuracy and the computation speed. In the following experiments, the parameters $\alpha$ and Q are set at $\alpha =0.5$ and $Q=3$.

Figure 4a,b verify the RMSE of DOA estimation versus SNR for white noise and colored noise, where $J=1000$. Consider three targets whose DOAs are ${\theta }_1=-8.2^{\circ }$, ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$. According to Figure 4a,b, we can conclude that the RSL0 algorithm achieves the best DOA estimation for white noise. For colored noise, the proposed algorithm and the RSL0 algorithm provide better angle estimation performance than the $l_1$-SVD, RV $l_1$-SVD and RV $l_1$-SRACV methods.

The processing time of signal reconstruction is compared in

Table 3. Average computation time of the signal reconstruction for different algorithms.

(M,N,P)	Average Computation Time (s)
	Proposed Method	${\mathit{l}}_1$-SVD	RV ${\mathit{l}}_1$-SVD	RV ${\mathit{l}}_1$-SRACV	RSL0
(4,4,2)	0.0269	2.1067	1.5499	1.3213	0.0135
(4,4,3)	0.0323	2.5139	2.0591	1.3368	0.0139
(5,5,2)	0.0273	2.4984	2.1071	2.1277	0.0252
(5,5,3)	0.0373	3.1170	2.7076	2.3842	0.0254
(6,6,2)	0.0328	2.6323	2.2911	3.2908	0.0438
(6,6,3)	0.0395	3.5439	3.1187	4.7397	0.0441

SVD (singular value decomposition), RV (real-valued), SRACV (sparse representation array covariance vectors), RSL0 (reweighted smoothed $l_0$-norm).

, where $J=1000$ and SNR = 0 dB. Different antenna numbers and target numbers are considered. For $p=2$, the DOAs are ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$. For $p=3$, they are ${\theta }_1=-8.2^{\circ }$, ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$. Table 3 verifies that, compared with the $l_1$-norm minimization based methods, the proposed algorithm has much lower computational complexity and is about two orders of magnitude faster.

Figure 5a,b illustrates the RMSE versus snapshots for Gaussian white noise and colored noise, where SNR = 0 dB, and the DOAs of three targets are ${\theta }_1=-8.2^{\circ }$, ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$. It can be observed that, for white noise, the RSL0 algorithm has the minimum RMSE. For colored noise, the proposed algorithm provides the best angle estimation performance when the snapshots vary from $J=450$ to $J=4050$, whereas, in this case, the RV $l_1$-SRACV method performs worse than the other algorithms.

Figure 6a,b verifies the DOA estimation performance of different sparse methods versus angle separation for Gaussian white noise and colored noise, where $J=1000$, and the DOAs of two targets are considered as ${\theta }_1=0^{\circ }$ and ${\theta }_2={\theta }_1+\tilde{\theta }$ with $\tilde{\theta }$ varying from $\tilde{\theta }=4^{\circ }$ to $\tilde{\theta }={16}^{\circ }$. For the white noise in Figure 6a, the RSL0 algorithm in [22] performs the best, and the proposed algorithm provides better estimation performance than the $l_1$-SVD, RV $l_1$-SVD and RV $l_1$-SRACV methods. For the colored noise in Figure 6b, the performance of the proposed algorithm stands out among the analyzed methods, while the performance of the RV $l_1$-SRACV is severely degraded.

Figure 7a,b verifies the relationships between the target resolution probability and SNR for different sparse methods, where $J=1000$, and three targets ${\theta }_1=-8.2^{\circ }$, ${\theta }_2=0^{\circ }$ and ${\theta }_3={21}^{\circ }$ are considered. The ith Monte Carlo trial is regarded as a successful one if the estimation results ${\overline{\theta }}_1$, ${\overline{\theta }}_2$ and ${\overline{\theta }}_3$ satisfy $\max \left\{\right|{\overline{\theta }}_1-{\theta }_1|,|{\overline{\theta }}_2-{\theta }_2|,|{\overline{\theta }}_3-{\theta }_3\left|\right\}\le 0.1^{\circ }$. Figure 7a,b verifies that the RSL0 algorithm performs the best for the white noise, and the proposed algorithm provides the highest target resolution probability for the colored noise. Colored noise causes the serious performance degradation for the analyzed methods except for the proposed algorithm.

6. Conclusions

Direction-of-arrival estimation is usually confronted with the MMV case. In this paper, for DOA estimation in MIMO radar, we have proposed a joint smoothed $l_0$-norm algorithm tailored for the multiple measurement vector case. By designing the new joint smoothed function and its gradient-based sparse signal reconstruction, the proposed method is a fast sparse DOA estimation algorithm that can solve the MMV problem. In addition, it is applicable to the colored Gaussian noise. The extended applications and the computational complexity of the proposed method are analyzed. The experimental results have verified that the proposed algorithm is about two orders of magnitude faster than the $l_1$-norm minimization based methods, such as $l_1$-SVD, RV $l_1$-SVD and RV $l_1$-SRACV, and performs well for both white and colored Gaussian noises.