 Research
 Open Access
 Published:
Twodimensional SLIM with application to pulse Doppler MIMO radars
EURASIP Journal on Advances in Signal Processing volume 2015, Article number: 69 (2015)
Abstract
A twodimensional (2D) sparse signal model is developed for pulse Doppler MIMO radars. Using this model, we develop the 2D sparse learning via iterative minimization (2D SLIM) algorithm. Simulation results show that the 2D SLIM compared to the 1D SLIM drastically reduces the computational burden while both of them have the same performance. Also, for estimation of rangeangleDoppler parameters, the 2D SLIM outperforms the matched filter (MF), smoothed L0norm (SL0), iterative adaptive approach (IAA), and spectral projected gradient for l _{1}norm minimization (SPGL1) algorithms.
Introduction
Multipleinput multipleoutput (MIMO) radars by exploiting multiple transmitters and receivers have recently been introduced [1–3]. It is well known that in this structure due to making use of orthogonal (or highly uncorrelated) transmit signals, the received signals can easily be separated. MIMO radars are often divided into two categories based on antenna placement. In the first one, the transmit and receive antennas are widely separated, and thus, the targets are observed from different directions dealing with target fluctuation fading [4–7]. In the second category, however, antennas are collocated so that the different phases from received signals can be extracted by the receivers. In this case, due to the waveforms diversity, a higher spatial resolution is achieved compared to the traditional radars. Also, in MIMO radars, target detection and parameter estimation are improved by suitably designing transmit beampattern [8–13]. Here, we consider the second structure.
From a sparsity perspective, in most radar applications, the number of targets located in the radar surveillance area is much smaller than the whole number of rangeangleDoppler bins. Thus, a sparse model can be derived for the received signal, and accordingly sparse signal recovery algorithms can be used for estimating the target parameters including range, Doppler frequency, and angle [14–18]. The aim of using sparse solution in a radar system is to more accurately estimate the target parameters compared to the traditional methods such as matched filters (MF).
Compressed sensing (CS), which is rooted on the principles of sparsity theory, has recently received considerable attention in MIMO radars [19–21]. Although, the main goal of the CS problem is to reduce the sampling rate lower than the Nyquist criterion, here we mainly focused on achieving accurate estimates for target parameters with much lower computations. For this purpose, an efficient technique is the sparse learning via iterative minimization (SLIM) algorithm which is computationally simple compared to the iterative adaptive approach (IAA) and focal underdetermined system solver (FOCUSS) algorithms due to the use of the conjugate gradient least squares (CGLS) algorithm [22]. This algorithm, which is a 1D algorithm (or namely 1D SLIM), is a maximum a posteriori (MAP) estimator which maximizes a posteriori Bayesian model. An important characteristic of 1D SLIM is incorporation of l _{ q }norm optimization (0 < q ≤ 1) in comparison with the l _{1}norm in order to reach sparser solutions and more accurate estimates. In [22], the 1D SLIM has been developed for MIMO radars by using only one pulse. Based on [22], this algorithm estimates a sparser vector compared to the l _{1}norm algorithm. On the other hand, the smoothed L0 (SL0) algorithm has been presented for twodimensional (2D) sparse problems [23]. In this algorithm, a discontinuous l _{0}norm function is approximated by a continuous one and then a sparse solution is reached using the steepest ascent algorithm followed by a projection onto a feasible set [24–27]. In [28], this algorithm has been applied to pulse Doppler radars with a lot of advantages such as target velocity extraction and pulse integration. However, this algorithm has been presented by using an approximated l _{0}norm function for which we will later show that it achieves a lower performance in sparse signal recovery at low signaltonoise ratios (SNRs) or for a small number of pulses compared to the SLIM algorithm. Also, the 2D IAA which is a nonparametric algorithm is presented in [29] for a general sparse solution. However, as it will be shown in the simulation section, its performance is poor at low SNRs and also for a small number of pulses.
In [30], a lowcomplexity CS approach is developed by decoupling the range, Doppler frequency, and angle parameters. It is assumed that the estimates of azimuth angles are obtained from one pulse by discretizing the angle space. Then, the Doppler estimates are extracted by combining the data of multiple pulses and using the initial estimated angles. Based on angleDoppler estimates, the range is then estimated using frequencyvarying received pulses. However, two problems are still of concern. First, the number of targets within the radar surveillance area is, in practice, so large that the angle space will not be sparse enough to apply the CS theory. In addition, a huge number of antennas are needed in a MIMO radar to discriminate among a large number of targets in the angle space with an acceptable resolution. Secondly, the SNR of one pulse is not sufficient for estimation of targets’ angles.
In this paper, we develop a 2D sparse model for pulse Doppler MIMO radar signals and find its relation with the Kronecker product factorization in the 1D model. To solve the 2D sparse signal equation, it can be converted to a 1D model and be recovered using 1D sparse recovery algorithms. However, this leads to a very large number of computations for which we will introduce here a new simpler technique. Therefore, a 2D SLIM algorithm is proposed for direct solution of a 2D sparse signal equation. In this approach, the Kronecker factorization is used to separate a largedimension matrix into two smaller matrices. This procedure leads to reducing the number of products and therefore decreasing the computation cost and required memories.
Moreover, we develop the 2D version of the wellknown 1D matched filter (1D MF) for comparison with the proposed 2D SLIM. Moreover, we compare the 2D SLIM algorithm with spectral projected gradient for l _{1}norm minimization (SPGL1) algorithm which is appropriate for largescale sparse recovery problems and complexvalued data [31].
The paper is organized as follows. In Section 2, in a common scenario for aircraft surveillance MIMO radar in which the intrapulse Doppler shift is negligible, a 2Dsparse model is developed. To solve this 2D sparse signal equation, in Section 3, the 2D SLIM and 2D MF algorithms are derived. The computational complexities of 1D and 2D SLIM algorithms are discussed in Section 4. Using simulation results in Section 5, computational complexities and the performances of different algorithms are compared, and Section 6 concludes the paper. The list of notations used in this work is shown in Table 1.
Signal model for pulse Doppler MIMO radars
Figure 1a shows a typical radio frequency (RF) transmit pulse train of a pulse Doppler radar in which τ is the pulse width. For the received signals, two different scenarios may be considered for extraction of the Doppler shift/frequency, f _{ d } [32, 33]. In a general scenario shown in Fig. 1b, τ and f _{ d } are large enough to have f _{ d } τ > 1, for which at least one period of f _{ d } lies within the receive pulse width. However, in the second scenario shown in Fig. 1c, we have f _{ d } τ < 1 (usually f _{ d } τ < < 1), in which case several pulses are required to extract f _{ d }. This scenario; which we have considered in this work, is commonly encountered in aircraft surveillance radars [32]. In this scenario, the effect of the Doppler frequency on each pulse is negligible and may be viewed as the sampling of the Doppler signal by the radar pulse repetition frequency (PRF), f _{ r }.
The data model and problem formulation for Fig. 1c are presented as follows. As shown in Fig. 2, the transmit signal is a train of N _{ P } probing pulses each of which containing N _{ s } subpulses with the bandwidth B. We assume that the targets are located behind the maximum unambiguous ranges with no ambiguity in the Doppler frequency interval of interest \( \left[\frac{{f}_r}{2},\frac{f_r}{2}\right) \). This interval is divided into N _{ D } Doppler bins as
Now, we calculate the amount of phase shift over one subpulse caused by the Doppler phenomenon. Thus, for the dth Doppler bin, the Doppler phase shift over one subpulse can defined as
By considering the waveform for subpulses \( {\boldsymbol{s}}_i\in {\mathrm{\mathbb{C}}}^{1\times {N}_s}, \) i = 1, …, M _{ t } as the code sequence of the ith transmit antenna, transmit signals are defined by an M _{ t } × N _{ s } matrix as
Also, the range dimension of surveillance area is divided into N _{ R } bins. Accordingly, the largest possible delay between the transmit and receive pulses is N _{ R } − 1. Then, the transmitted pulse waveforms can be arranged into the matrix \( \tilde{\boldsymbol{S}}, \) so that
where \( {\mathbf{0}}_{M_t\times \left({N}_R1\right)} \) is an M _{ t } × (N _{ R } − 1) matrix of zeros, and we have \( \tilde{\boldsymbol{S}}\in {\mathrm{\mathbb{C}}}^{M_t\times \left({N}_s+{N}_R1\right)} \). Also, we assume that the angular interval of interest θ _{ a } is divided into N _{ A } angular bins (a = 1, ⋯, N _{ A } ). Then, in a uniform linear array, the steering vectors of M _{ t } transmit and M _{ r } receive antennas are respectively denoted by
and
where λ _{0} is the radar carrier wavelength and Δ _{ t } and Δ _{ r } show the distances between two adjacent transmit and receive antennas, respectively. Therefore, the pth received pulse matrix \( {\boldsymbol{Y}}_p\in {\mathrm{\mathbb{C}}}^{M_r\times \left({N}_s+{N}_R1\right)} \) can be written as
where Ε _{ p } is the additive white Gaussian noise matrix for the pth pulse (p = 1, 2, ⋯, N _{ P }), T _{ R } is the ratio of pulse repetition interval (PRI) over the single subpulse duration, and
is an (N _{ s } + N _{ R } − 1) × (N _{ s } + N _{ R } − 1) shift matrix for the rth range bin. Also, α _{ r,a,d } for r = 1, ⋯, N _{ R }, a = 1, ⋯, N _{ A }, and d = 1, ⋯, N _{ D } denote the return complex reflection coefficients of targets corresponding to the radar crosssection. In this scenario, since the number of rangeangleDoppler bins, in which actual targets are detected, is much smaller than the total number of radar bins, most of reflection coefficients are zero, and thus, the radar signal model can be assumed sparse. We assume that complex reflection coefficients during pulse repetitions are constant.
We will show that (7) can be converted to a 2D sparse from in which the unknown parameters α _{ r,a,d } are presented in matrix form. By defining and using
and
in (7), we obtain
where e _{ p } = vec(E _{ p }) is a complex Gaussian noise vector with zero mean and covariance matrix I. Equivalently, in a more compact form, we get
where \( \boldsymbol{X}=\left[\begin{array}{ccc}\hfill {\boldsymbol{x}}_1\hfill & \hfill \cdots \hfill & \hfill {\boldsymbol{x}}_{N_D}\hfill \end{array}\right] \) contains the complex reflection coefficients corresponding to the radar crosssection, and \( {\boldsymbol{\theta}}_p={\left[\begin{array}{ccc}\hfill {e}^{j{T}_R\left(p1\right){\omega}_1}\hfill & \hfill \cdots \hfill & \hfill {e}^{j{T}_R\left(p1\right){\omega}_{N_D}}\hfill \end{array}\right]}^T. \)
Next, by defining y _{ p }, p = 1, 2, ⋯, N _{ P } as the columns of matrix Y, we have
where \( \boldsymbol{E}=\left[\begin{array}{ccc}\hfill {\boldsymbol{e}}_1\hfill & \hfill \cdots \hfill & \hfill {\boldsymbol{e}}_{N_P}\hfill \end{array}\right] \) and \( \boldsymbol{\Theta} =\left[\begin{array}{ccc}\hfill {\boldsymbol{\theta}}_1\hfill & \hfill \cdots \hfill & \hfill {\boldsymbol{\theta}}_{N_P}\hfill \end{array}\right] \).
Equation (14) presents a 2D sparse signal model for pulse Doppler MIMO radars where \( \boldsymbol{Y}\in {\mathrm{\mathbb{C}}}^{M_r\left({N}_s+{N}_R1\right)\times {N}_P}, \) \( \boldsymbol{A}\in {\mathrm{\mathbb{C}}}^{\left[{M}_r\left({N}_s+{N}_R1\right)\right]\times \left[{N}_R{N}_A\right]}, \) \( \boldsymbol{\Theta} \in {\mathrm{\mathbb{C}}}^{N_D\times {N}_P}, \) and \( \boldsymbol{X}\in {\mathrm{\mathbb{C}}}^{N_R{N}_A\times {N}_D} \). Due to the underdetermined nature of (14) for a sparse model (i.e., M _{ r }(N _{ s } + N _{ R } − 1) < N _{ R } N _{ A } and N _{ P } < N _{ D }), it has no unique solution. Our goal is to find the sparsest matrix for X in which we have as many zero elements as possible.
The 2D sparse signal model given by (14) can be converted to a 1D model by using the following property [34]:
Therefore, we have
where x = vec(X), y = vec(Y), e = vec(E), and Φ = Θ ^{T} ⊗ A. Although x can be computed using 1D sparse recovery algorithms such as IAA [35] and SLIM [22], due to the large dimension of \( \boldsymbol{\Phi} \in {\mathrm{\mathbb{C}}}^{\left[{N}_P{M}_r\left({N}_s+{N}_R1\right)\right]\times \left[{N}_R{N}_A{N}_D\right]}, \) 1D solutions are computationally extremely expensive. Accordingly, due to the smaller number of products appeared in (14) compared to (16), developing the 2D algorithm for direct solution of (14) leads to an extreme reduction of computational load compared to 1D one.
2D sparse signal recovery
In this section, at first, we give an overview of 1D SLIM algorithm. Then, a 2D SLIM algorithm is proposed for pulse Doppler MIMO radars by direct solution of (14), which leads to a lower computational cost compared to the 1D algorithms. In addition, for comparison purposes, the 2D SLIM is compared with the 2D SL0 [23] and 2D IAA [29] algorithms recently introduced for 2D sparse recovery problems. Moreover, we develop the 2D MF for comparison purposes.
Overview on 1D SLIM algorithm
In the 1D SLIM algorithm which is based on a MAP approach, the unknown complex reflection coefficients are considered as random variables with the following prior distribution
where q is a free parameter between 0 < q ≤ 1, x _{ n } is the nth component of vector x, and N = N _{ R } N _{ A } N _{ D }. The above prior can lead to a more accurate estimation of the sparse vector x. When q → 0, this prior will have a sharper peak at 0 and thus causes sparser estimates in the Bayesian inference. For q = 1, this prior will be similar to the Laplace prior \( f\left(\boldsymbol{x}\right)\propto {e}^{{\left\Vert \boldsymbol{x}\right\Vert}_1} \) with finite peaks at 0. By assuming an independent and complex Gaussian distribution with zero mean and variance η I for additive noise e, the conditional distribution for measurement data vector y can be defined as \( f\left(\boldsymbol{y}\Big\boldsymbol{x},\eta \right)=\mathcal{C}\mathcal{N}\left(\boldsymbol{Ax},\ \eta \boldsymbol{I}\right) \). By assuming a uniform prior for η as f(η) ∝ 1, the cost function for the SLIM based on the MAP approach is defined as
where K = M _{ r } N _{ P }(N _{ s } + N _{ R } − 1). By taking the negative of logarithm of (18), the cost function will be equivalent to
By minimizing (19) with respect to x and η and using a heuristic approach, an iterative solution is obtained for the 1D SLIM in two steps [22]:

1.
Iterative estimation of sparse vector x,
where the superscript ^{t} shows the iteration number, Π = diag{ϑ}, and ϑ _{ n } = x _{ n }^{2 − q}, n = 1, ⋯, N _{ R } N _{ A } N _{ D }, ϑ _{ n } is the nth component of vector ϑ.

2.
Iterative estimation of noise power η,
2D SLIM
Our aim in this subsection is to develop the 2D version of (20) and (21) using Φ = Θ ^{T} ⊗ A and (14) and following Hadamard matrix property:
Since matrix inversion in (20) has an extreme computational load, we propose to use the conjugate gradient least square (CGLS) algorithm [36] to solve (14) with lower computations. To do so, by defining the vector u = (ΦΠΦ ^{H} + η I)^{− 1} y, the 2D CGLS is derived in the Appendix as summarized in Table 2. Note that the superscript ^{t} has been suppressed for simplicity and the components of matrix Σ are Σ_{ ij } = Γ_{ ij } ^{1/2} where
To derive the 2D SLIM, the relationship between Π and Γ can be presented as
where the elements of Γ are given by (23). By substituting (24) and Φ = Θ ^{T} ⊗ A in (20), we obtain
Then, by using the Kronecker product property [37]
and using (15) and (22), (25) is shown as
where \( \boldsymbol{U}\in {\mathrm{\mathbb{C}}}^{M_r\left({N}_s+{N}_R1\right)\times {N}_P} \) is obtained by the 2D CGLS algorithm as presented in Table 2 and u = vec(U). Note that X = Γ ⊙ (A ^{H} U Θ ^{H}), and thus for the (t + 1)th iteration of step 1 of the 2D SLIM, we have
Also, from (14) and (16), we can easily express step 2 of the 2D SLIM as
For initialization of the 2D SLIM, i.e., X ^{(0)}, we make use of minimum l _{2}norm solution using the economysize QR decomposition [28] by incorporation of a hard threshold operator in order to sparsify X ^{(0)}. The required criterion for stopping the 2D SLIM recursion is developed based on the 1D version [22] as
where Δ is a small positive constant.
The parameter q plays an important role in sparse signals recovery using the SLIM. For q → 0, it achieves a sparser solution compared to the case when q → 1; however, it is hard to determine the sparsity level. To cope with this problem, q can be estimated appropriately by the Bayesian information criterion (BIC) [22] given by
where h(q) is the number of selected peaks in the output of the SLIM algorithm that is executed for a particular value of q. To explain how we choose the number of selected peaks (h(q)), at first, the absolute value of the SLIM algorithm output defined as X is sorted in a descending order. Then, the largest peak is selected and the other values of matrix X are set to zero to form the matrix \( \widehat{\boldsymbol{X}} \). Using (28), the BIC is computed for \( \widehat{\boldsymbol{X}} \) and h(q) = 1. In the round, the two largest peaks are selected and the other values of X are set to zero, and the related BIC is computed for \( \widehat{\boldsymbol{X}} \) and h(q) = 2 and so on. The value of h(q) is the number of the selected peaks that yields the lowest BIC. After running the 2D SLIM for a selected set of q and computing h(q), we choose that q which minimizes the BIC. The factor 5 in (28) shows the number of unknown parameters including range, angle, Doppler frequency, and the complex reflection coefficients of targets.
2D MF
The 1D MF output has already been derived for MIMO radars as [22]
Noting that the columns of Φ are defined by the vectors \( {\left(\boldsymbol{\Phi} \right)}_n^c,\ n=1,2,\dots, {N}_R{N}_A{N}_D \) and using Φ = Θ ^{T} ⊗ A and (26), the numerator of (32) is obtained as
Next, the denominator of (32) is shown by using the Kronecker product properties as
where n = (j − 1)N _{ R } N _{ A } + i, i = 1, ⋯, N _{ R } N _{ A }, and j = 1, ⋯, N _{ D }. By substituting (33) and (34) in (32) and converting to the 2D form, we obtain the 2D MF as
where the elements of \( \boldsymbol{\Lambda} \in {\mathrm{\mathbb{C}}}^{N_R{N}_A\times {N}_D} \) are defined by Λ_{ ij } and ⊘ is the elementwise matrix division.
Computational complexity of 1D and 2D SLIM
In the 1D SLIM, the main computational cost in each iteration belongs to the product of Φ and a vector like x. For the 2D SLIM, this product is converted to a 2D form as AX Θ using the Kronecker factorization Φ = Θ ^{T} ⊗ A.
The main difference between the 1D and 2D forms from a computational point of view is in the number of flops for calculating Φ x and AX Θ. The complexities of Φ x and AX Θ are O(N _{ P } M _{ r }(N _{ s } + N _{ R } − 1)N _{ R } N _{ A } N _{ D }) and O(M _{ r }(N _{ s } + N _{ R } − 1)N _{ R } N _{ A } N _{ D }) + O(N _{ P } M _{ r }(N _{ s } + N _{ R } − 1)N _{ D }), respectively. We have assumed that the product matrix AX is computed first, and then the result is multiplied by Θ. By comparing these two computational complexities, it is demonstrated that the ratio of the 2D processing load over that of its equivalent 1D is \( \frac{1}{N_P}+\frac{1}{N_R{N}_A} \). If it is assumed N _{ P } ≪ N _{ R } N _{ A }, then this ratio is equal to \( \frac{1}{N_P} \).
In addition, the Kronecker factorization can take advantage of multicore processors [38] by which the 2D SLIM algorithm can be parallelized and solved. Then, the speed of this algorithm compared to the 1D one is approximately increased by a factor proportional to the number of processing cores.
Numerical examples
To show the computational efficiency of the 2D SLIM, which is the main goal of this work, we compare the running time of the 1D version of SLIM, IAA, SL0, and MF algorithms with their 2D versions. To make a fair comparison between 1D and 2D SLIM algorithms, we use the CGLS in the 1D SLIM to reduce the computational cost of matrix inversion. Moreover, we express the product of the diagonal matrix Π and the vector Φ ^{H} u by the Hadamard product as x = ΠΦ ^{H} u = ϑ ⊙ (Φ ^{H} u) and do a similar procedure for the CGLS steps.
Also, the performance of the above algorithms for estimation of rangeangleDoppler parameters in pulse Doppler MIMO radars is compared to that of the SPGL1, which is an l _{1}based approach. Note that since the performance analysis of 1D and 2D algorithms is similar, in the related figures, we will ignore repeating 1D or 2D terms. The noise vectors e _{ p }, p = 1, ⋯, N _{ P } are mutually independent and each vector consists of zero mean complex white Gaussian noise with covariance matrix σ ^{2} I. The signaltonoise ratio (SNR) is defined for each target located at the (r,a,d)th rangeangleDoppler bin as
The number of targets is N _{ t } = 40 and the SNR of each target is 10 dB. We consider the transmit signal with a cyclic approach [39] with N _{ s } = 32. The transmit and receive antennas are uniform linear arrays with Δ_{ t } = 2.5λ _{0}, Δ_{ r } = 0.5λ _{0}, and M _{ t } = M _{ r } = 5. The carrier frequency, subpulses bandwidth, and the PRF are respectively f _{ c } = 1 GHz, B = 10 MHz, and f _{ r } = 2 kHz. The surveillance field is divided into N _{ R } = 20 range bins, N _{ A } = 31 angular bins between − 30° to 30° with 2° angular resolution, and N _{ D } = 40 Doppler bins with the resolution of 50 Hz. The lower and upper thresholds for limiting the amplitude of output signals are − 30 and 10 dB, respectively.
The SL0, IAA, and SLIM algorithms are initialized based on an l _{2}norm solution with incorporating a hard threshold set to − 20 dB. It means that when the absolute value of each component of initial X ^{(0)} is lower than − 20 dB, it is set to 0. The free parameters for 1D and 2D SL0 are μ = 1, ε = 0.002, σ _{0} = 0.05, ρ = 0.5, and K = 10, and for the 1D and 2D SLIM, we have Δ = 0.01. Also, the stopping parameter for the CGLS is set to ε = 0.05. BIC values versus q are shown in Fig. 3 for N _{ P } = 5. As seen, q = 0.1 can be chosen for this scenario.
We consider two different scenarios in which the targets may fall onto or off the grid points. In Figs. 4 and 5 where the targets fall onto the grid points, the targets’ true locations are shown with circles and colorcoded rectangles which correspond to the estimated amplitudes in dB. As seen in Figs. 4a and 5a, the MF approach suffers from accurate estimation of the targets’ anglerangeDoppler parameters due to appearing large side lobe levels. Also, in Figs. 4b and 5b, the SL0 does not perform acceptable results when the number of pulses is small. Instead, the SLIM, IAA, and SPGL1 have mostly estimated the range, angle, Doppler frequency parameters; however, the former algorithm is more accurate. This is due to the fact that the fewer the number of pulses are, the more underdetermined (16) will be. To cope with such cases, then, a more effective sparse recovery algorithm is required.
The second scenario, in which the targets do not fall onto the grid points, is more practical. In this scenario, we choose the targets randomly from all possible bins and for each target, a random number proportional to the angle resolution and a random number proportional to the Doppler resolution are added to the angle and Doppler of the target. In this case, the targets will fall off the grid points. In this simulation, we consider N _{ t } = 20 and SNR = 0 dB. Figure 6 demonstrates that the SLIM, IAA, and SPGL1 have mostly captured the targets that fall out of grid points, while the MF and SL0 have failed. Also, the SLIM algorithm has better performance compared to other ones.
In the next experiment, using Monte Carlo simulations, we compare the mean squared error (MSE) and peaktoripple ratio (PRR) of different algorithms. The MSE of target scene recovery is defined as \( MSE={\left\Vert \widehat{X}\boldsymbol{X}\right\Vert}_F^2/N, \) where N = N _{ R } N _{ A } N _{ D }. Also, for N _{ t } targets located at the angleDoppler range bins {(k _{ i }, l _{ i }), i = 1, …, N _{ t }}, the PRR is given by
The results are presented based on 1000 independent trials. The MSEs and PRRs are respectively shown versus the number of pulses in Fig. 7a, d for N _{ t } = 50 and SNR = 10 dB, versus the number of targets in Fig. 7b, e for N _{ P } = 5 and SNR = 10 dB, and versus the SNR in Fig. 7c, f for N _{ P } = 8 and N _{ t } = 50. By comparing the above results, one can clearly see that the SLIM algorithm achieves a lower MSE and a larger PRR in estimation of rangeangleDoppler parameters.
Also, we calculate the receiver operating characteristic (ROC) for different algorithms in order to compare their detection performance. To obtain the ROC curve, we select a threshold τ _{ i }. Any local maximum of the absolute value of X that is larger than τ _{ i } will be considered as a target. The threshold τ _{ i } is then varied within an interval [τ _{ L } τ _{ H }] and for each τ _{ i }, the number of detected actual or false targets is recorded. By repeating 1000 trials of the experiment, we calculate the probability of detection P _{ d } of actual targets, and the probability of false alarm P _{ f } of false targets for different values of τ _{ i } as
In Fig. 8, we compare the ROC of different algorithms for N _{ P } = 8, N _{ t } = 200, and SNR = 10 dB. As seen, in all cases the probability of detection of the SLIM is better than that of the other ones.
The running times of different algorithms are presented in Table 3 for N _{ P } = 20 by averaging over 100 independent trials of the experiments. We run our MATLAB 8.1 files on a PC with an Intel Core i5/3.4 GHz with 4 GB memory.
Obviously, in all cases, the 2D versions are much more efficient than 1D ones. Correspondingly, in Fig. 9, the running times are depicted as a function of N _{ P }. One can see that by increasing the number of pulses, the running times of 2D algorithms remain almost constant. In other words, the computational costs of the latter algorithms are less sensitive to the number of pulses, N _{ P }. The reason is that the computational complexity of 2D SLIM does not depend on N _{ P } as explained in Section 4. As a result, according to Fig. 7a, d, in order to increase the performance of 2D algorithms in estimation of radar parameters, we can easily use a larger N _{ P }, while for 1D algorithms, this leads to a large increase in the running time. Moreover, it is seen from Fig. 9 that the running times of 1D algorithms are much longer than those of the 2D ones, revealing their higher computational complexities. Note that although the 2D MF has the least running time, it also generates the lowest detection and recovery performance as already depicted in Figs. 7 and 8.
Conclusions
We derived the 1D and 2D sparse signal models for pulse Doppler MIMO radars. The 2D SLIM algorithm was derived by direction solution of the 2D sparse model. Due to using a lower number of products in the corresponding relationships, the computational cost of 2D SLIM compared to that of the 1D one was extremely reduced. Also, simulation results showed that the 2D SLIM outperforms the other algorithms in accurate estimation of range, angle, and Doppler parameters.
References
 1.
E Fishler, A Haimovich, R Blum, D Chizhik, L Cimini, R Valenzuela, MIMO Radar: An Idea Whose Time has Come. Proceedings of the IEEE Radar Conference, 2004, pp. 71–78
 2.
FC Robey, S Coutts, D Weikle, JC McHarg, and K Cuomo, in MIMO Radar Theory and Experimental Results. Proceedings of 38th Asilomar Conference on Signals, Systems and Computers, vol 1 (IEEE, Pacific Grove, CA, 2004), pp. 300–304
 3.
J Li, P Stoica, MIMO radar—diversity means superiority. Proc. 4th Adaptive Sensor Array Process Workshop. (2006). doi:10.1002/9780470391488
 4.
N Lehmann, F Fishler, AM Haimovich, RS Blum, D Chizhik, L Cimin, et al. Evaluation of Transmit Diversity in MIMORadar Direction Finding. IEEE Transactions on Signal Processing, 2007, pp. 2215–2225
 5.
E Fishler, A Haimovich, R Blum, L Cimini, D Chizhik, and R Valenzuela, in Performance of MIMO Radar Systems: Advantages of Angular Diversity. Proceedings of 38th Asilomar Conference on Signals, Systems and Computers, vol 1 (IEEE, Pacific Grove, CA, 2004), pp. 305–309
 6.
AD Maio, M Lops, Design principles of MIMO radar detectors. IEEE Trans. Aerosp. Electron. Syst. 43(3), 886–898 (2007)
 7.
A Sheikhi, A Zamani, Temporal coherent adaptive target detection for multiinput multioutput radars in clutter. IET Trans. Radar Sonar Navig. 2(2), 86–96 (2008)
 8.
J Li, P Stoica, MIMO radar with collocated antennas: review of some recent work. IEEE Signal Process. Mag. 24(5), 106–114 (2007)
 9.
D W Bliss, KW Forsythe, in MultipleInput MultipleOutput (MIMO) Radar and Imaging: Degrees of Freedom and Resolution. 37th Asilomar Conference on Signals, Systems and Computers, vol 1 (IEEE, Pacific Grove, CA, 2003), pp. 54–59
 10.
K Forsythe, D Bliss, G Fawcett, in MultipleInput MultipleOutput (MIMO) Radar: Performance Issues. 38th Asilomar Conference on Signals, Systems and Computers, vol 1, (IEEE, Pacific Grove, CA, 2004), pp. 310–315
 11.
J Li, P Stoica, MIMO radar—diversity means superiority. The Fourteenth Annual Workshop on Adaptive Sensor Array Processing (invited). (2006). doi:10.1002/9780470391488.
 12.
P Stoica, J Li, Y Xie, On probing signal design for MIMO radar. IEEE Trans. Signal Process. 55(8), 4151–4161 (2007)
 13.
X Luzhou, J Li, P Stoica, Target detection and parameter estimation for MIMO radar systems. IEEE Trans. Aerosp. Electron. Syst. 44(3), 927–939 (2008)
 14.
RG Baraniuk, Compressive Sensing. IEEE Signal Processing Magazine, 2007, pp. 118–121
 15.
DL Donoho, Compressed Sensing. IEEE Transactions on Information Theory, 2006, pp. 1289–1306
 16.
EJ Candès, Compressive Sampling, in Proc. International Congress of Mathematicians, 2006, pp. 1433–1452
 17.
EJ Candès, MB Wakin, An introduction to compressive sampling. IEEE Signal Process. Mag. 25(2), 21–30 (2008)
 18.
MA Herman, T Strohmer, Highresolution radar via compressed sensing. IEEE Trans. Signal Proc. 57(6), 2275–2284 (2009)
 19.
AP Petropulu, Y Yu, HV Poor, Distributed MIMO radar using compressive sampling. in Proc, 42nd Asilomar Conf, Signals, Syst, Comput., (IEEE, Pacific Grove, CA, 2008), pp. 41–44
 20.
T Strohmer, B Friedlander, Compressed sensing for MIMO radar—algorithms and performance. in Proc. 43rd Asilomar Conf. Signals, Syst. Comput., (IEEE, Pacific Grove, CA, 2009), pp. 464–468
 21.
Y Yu, AP Petropulu, HV Poor, MIMO radar using compressive sampling. IEEE J. Sel. Topics Signal Process. 4(1), 146–163 (2010)
 22.
X Tan, W Roberts, J Li, P Stoica, Sparse learning via iterative minimization with application to MIMO radar imaging. IEEE Trans. Signal Process. 59(3), 1088–1101 (2011)
 23.
A Ghaffari, M BabaieZadeh, C Jutten, Sparse decomposition of two dimensional signals. in IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP, (IEEE, Taipei, 2009), pp. 3157–3160
 24.
A Eftekhari, M BabaieZadeh, C Jutten, H Moghaddam, RobustSL0 for stable sparse representation in noisy settings. in IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP, (IEEE, Taipei, 2009), pp. 3433–3436
 25.
H Mohimani, M BabaieZadeh, C Jutten, A Fast Approach for Overcomplete Sparse Decomposition Based on Smoothed l _{ 0 } Norm. IEEE Transactions on Signal Processing, 2009, pp. 289–301
 26.
GH Mohammadi, M BabaieZadeh, and C Jutten, Complexvalued sparse representation based on smoothed l0 norm. in IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP, (IEEE, Las Vegas, 2008), pp. 3881–3884
 27.
M Hyder, K Mahata, An improved smoothed l _{ 0 } approximation algorithm for sparse representation. IEEE Trans. Signal Process. 58(4), 2194–2205 (2010)
 28.
M Hyder, K Mahata, RangeDoppler Imaging via Sparse Representation. Radar Conference, 2011, pp. 486–491
 29.
M JabbarianJahromi, MH Kahaei, Twodimensional iterative adaptive approach for a sparse matrix solution. Electron. Lett. 50, 45–47 (2014)
 30.
Y Yu, AP Petropulu, HV Poor, Reduced Complexity AngleDopplerRange Estimation for MIMO Radar That Employs Compressive Sensing. Proc. 43nd Asilomar Conference on Signals, Systems and Computers, 2009
 31.
E Van Den Berg, MP Friedlander, SPGL1: a solver for largescale sparse reconstruction. http://www.math.ucdavis.edu/~mpf/spgl1/. Accessed June 2007.
 32.
M Skolnik, Introduction to Radar System, 3rd edn. (McGrawHill, New York, 2002)
 33.
M Skolnik, Radar Handbook, 3rd edn. (McGrawHill, New York, 2008)
 34.
KB Petersen, MS Pedersen, The Matrix Cookbook, 2008. Available: http://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf
 35.
T Yardibi, J Li, P Stoica, M Xue, AB Baggeroer, Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares. IEEE Trans. Aerosp. Electron. Syst. 46, 425–443 (2010)
 36.
R Aster, B Borchers, C Thurber, Parameter Estimation and Inverse Problem (Elsevier Academic, Burlington, MA, 2005)
 37.
AN Langville, WJ Stewart, The Kronecker Product and Stochastic Automata Networks. Elsevier Journal of computational and Applied Mathematics, 2004, pp. 429–447
 38.
G Wu, Z Zhang, E Chang, Kronecker Factorization for Speeding up Kernel Machines. SIAM International Conference on Data Mining, 2005, pp. 611–615
 39.
H He, P Stoica, J Li, Designing unimodular sequence sets with good correlations—including an application to MIMO radar. IEEE Trans. Signal Process. 57, 4391–4405 (2009)
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
About the authors
Mohammad Jabbarian Jahromi received the BS degree in electrical engineering from Shiraz University, Shiraz, Iran, in 2005 and the MSc degree from the Electrical and Electronic Engineering University Complex, Tehran, Iran, in 2008. He is currently working toward the PhD degree in the School of Electrical Engineering, Iran University of Science and Technology. From 2009, he was a research staff in the Signal & System Modeling laboratory, IUST. His current research interests are in the areas of array signal processing, radar signal processing, and sparse reconstruction.
Mohammad Hossein Kahaei received the BSc degree from Isfahan University of Technology, Isfahan, Iran, in 1986, the MSc degree from the University of the Ryukyus, Okinawa, Japan, in 1994, and the PhD degree in signal processing from the School of Electrical and Electronic Systems Engineering, Queensland University of Technology, Brisbane, Australia, in 1998. Since 1999, he has been with the School of Electrical Engineering, Iran University of Science and Technology, Tehran, Iran, where he is currently an Associate Professor and the head of Signal and System Modeling laboratory. His research interests include array signal processing with primary emphasis on compressed sensing, blind source separation, localization, tracking, and DOA estimation, and wireless sensor networks.
Appendix
Appendix
Developing 2d CGLS
We can convert the inversing term appeared in u to a least squares (LS) form as
where \( \boldsymbol{B}=\left(\begin{array}{c}\hfill {\boldsymbol{\Pi}}^{1/2}{\boldsymbol{\Phi}}^H\hfill \\ {}\hfill {\eta}^{1/2}\boldsymbol{I}\hfill \end{array}\right) \) and \( \boldsymbol{z}=\left(\begin{array}{c}\hfill 0\hfill \\ {}\hfill {\eta}^{1/2}\boldsymbol{y}\hfill \end{array}\right) \).
The solution of u is equivalent to minimization of the LS cost function \( {\left\Vert \boldsymbol{B}\boldsymbol{u}\boldsymbol{z}\right\Vert}_2^2 \) with respect to u, and it can be solved using the CGLS algorithm with the following steps [36]:

1.
\( {\alpha}_l={\boldsymbol{r}}_l^H{\boldsymbol{r}}_l/\left(\left({\boldsymbol{p}}_l^H{\boldsymbol{B}}^H\right)\left(\boldsymbol{B}{\boldsymbol{p}}_l\right)\right) \)

2.
u _{ l + 1} = u _{ l } + α _{ l } p _{ l }

3.
s _{ l + 1} = s _{ l } + α _{ l } Bp _{ l }

4.
r _{ l + 1} = B ^{H} s _{ l + 1}

5.
\( {\beta}_l={\boldsymbol{r}}_{l+1}^H{\boldsymbol{r}}_{l+1}/\left({\boldsymbol{r}}_l^H{\boldsymbol{r}}_l\right) \)

6.
p _{ l + 1} = − r _{ l + 1} + β _{ l } p _{ l }
where initializations are u _{0} = 0, s _{0} = − z, r _{0} = B ^{H} s _{0}, and p _{0} = − r _{0}. The stopping condition for CGLS is \( \frac{1}{N_z}{\left\Vert \boldsymbol{B}\boldsymbol{u}\boldsymbol{z}\right\Vert}_2<\varepsilon \), where N _{ z } is the number of elements of z.
To derive the 2D CGLS, we first present Bp _{ l } as
where vec(W _{ l }) = Π ^{1/2} Φ ^{H} p _{ l } with Π ^{1/2} = diag{vec(Σ)} and the elements of Σ are defined in (23) as Σ_{ ij } = Γ_{ ij } ^{1/2}. By using Φ = Θ ^{T} ⊗ A in vec(W _{ l }), we obtain
By incorporation of (15) and (22) in (42), we get vec(W _{ l }) = vec(Γ ⊙ (A ^{H} P _{ l } Θ ^{H})) where
Furthermore, V _{ l } in (41) is presented by V _{ l } = η ^{1/2} P _{ l } where p _{ l } = vec(P _{ l }).
Next, the term B ^{H} s _{ l + 1} of the CGLS is demonstrated using the Kroneker and Hadamard products properties shown in (15) and (22), and Φ = Θ ^{T} ⊗ A as
where \( {\boldsymbol{s}}_l=\left(\begin{array}{c}\hfill {\boldsymbol{g}}_l\hfill \\ {}\hfill {\boldsymbol{t}}_l\hfill \end{array}\right) \) with g _{ l } = vec(G _{ l }) and t _{ l } = vec(T _{ l }). By substituting the above results in the third and fourth steps of the CGLS and also presenting the vectors u _{ l } and r _{ l } in 2D forms as u _{ l } = vec(U _{ l }) and r _{ l } = vec(R _{ l }), respectively, the 2D CGLS is derived as summarized in Table 2. After some mathematical manipulations using the Hadamard and Kronecker product properties, the stopping condition for the 2D CGLS is obtained by extending \( \frac{1}{N_z}{\left\Vert \boldsymbol{B}\boldsymbol{u}\mathrm{z}\right\Vert}_2<\varepsilon \) as
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0), which permits use, duplication, adaptation, distribution, and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
JabbarianJahromi, M., Kahaei, M.H. Twodimensional SLIM with application to pulse Doppler MIMO radars. EURASIP J. Adv. Signal Process. 2015, 69 (2015). https://doi.org/10.1186/s1363401502546
Received:
Accepted:
Published:
Keywords
 Pulse Doppler MIMO radar
 Sparse learning via iterative minimization
 Twodimensional sparse signal model