Splitting Matching Pursuit Method for Reconstructing Sparse Signal in Compressed Sensing

Abstract

In this paper, a novel method named as splitting matching pursuit (SMP) is proposed to reconstruct 𝐾-sparse signal in compressed sensing. The proposed method selects 𝐹𝑙  (𝐹𝑙>2𝐾) largest components of the correlation vector 𝑐, which are divided into 𝐹 split sets with equal length 𝑙. The searching area is thus expanded to incorporate more candidate components, which increases the probability of finding the true components at one iteration. The proposed method does not require the sparsity level 𝐾 to be known in prior. The Merging, Estimation and Pruning steps are carried out for each split set independently, which makes it especially suitable for parallel computation. The proposed SMP method is then extended to more practical condition, e.g. the direction of arrival (DOA) estimation problem in phased array radar system using compressed sensing. Numerical simulations show that the proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods. The proposed method also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon, APES (amplitude and phase estimation) and GLRT (generalized likelihood ratio test) based on hundreds of snapshots.

1. Introduction

The standard noiseless model in compressed sensing is 𝑦=Φ𝑥, (1) where 𝑥∈ℜ^𝑁 is a 𝐾-sparse signal (𝐾≪𝑁), 𝑦∈ℜ^𝑀 is a measurement of 𝑥, and Φ is an 𝑀×𝑁 sensing matrix. The compressed sensing recovery problem is defined as follows: given 𝑦 and Φ, find a signal 𝑥 within the class of interest satisfies (1) exactly. The compressed sensing recovery process consists of a search for the sparsest signal 𝑥 that yields the measurement 𝑦. By defining the 𝑙₀ “norm” of a vector ‖𝑤‖₀ as the number of nonzero entries in 𝑤, the simplest way to pose a recovery algorithm is using the optimization𝑥=argmin_𝑦=Φ𝑤‖𝑤‖₀. (2)

Solutions to (2) therefore lead to algorithms for recovering 𝐾-sparse signals from 𝑀 linear measurements. In general, the minimization of (2) is NP-hard. An alternative to the 𝑙₀ “norm” used in (2) is to use the 𝑙₁ “norm”, defined as ‖𝑤‖₁=∑^𝑁_𝑛=1|𝑤(𝑛)|. The resulting adaptation of (2), known as basis pursuit (BP) [1], is formally defined as𝑥=argmin_𝑦=Φ𝑤‖𝑤‖₁. (3)

Since the 𝑙₁ “norm” is convex, (3) can be seen as a convex relaxation of (2). The optimization (3) can be modified to allow for noise in the measurements 𝑦=Φ𝑥+𝑒, where 𝑒 denotes an 𝑀×1 measurement noise vector. We simply change the constraint on the solution to 𝑥=argmin_{𝑦−Φ𝑤2≤𝜖}‖𝑤‖₁, (4) where 𝜖>‖𝑒‖₂ is an appropriately chosen bound on the noise magnitude. This modified optimization is known as basis pursuit with inequality constraints (BPIC) and is a quadratic program with polynomial complexity solvers [2]. The Lagrangian relaxation of this quadratic program is written as 𝑥=argmin‖𝑤‖₁+𝜆𝑦−Φ𝑤₂ (5) and is known as basis pursuit denoising (BPDN). There exist many efficient solvers to find BP, BPIC, and BPDN solutions; for an overview, see [3]. Unfortunately, the complexity of the linear programming algorithms for solving (3)~(5) is highly impractical for large-scale applications.

An alternative approach to sparse signal recovery is based on the idea of iterative greedy pursuit. The basic greedy algorithm is the matching pursuit (MP) [4]. OMP [5] is a variation of MP method, which adds a least-squares minimization step to MP method to obtain the best approximation over the chosen atoms. Unlike MP and OMP choosing just one atom at each time, ROMP, StOMP, SP, CoSaMP, and BAOMP methods choose several atoms at each iteration. Furthermore, regularized OMP (ROMP) [6], subspace pursuit (SP) [7], compressive sampling matching pursuit (CoSaMP) [8], and backtracking-based matching pursuit (BAOMP) [9] methods also use a two-step selection technique to carefully choose the atoms. While SP and CoSaMP have offered comparable theoretical reconstruction quality to the linear programming methods along with low reconstruction complexity, they require the sparsity level to be known for exact recovery. As an improvement, the BAOMP algorithm achieves the blind sparse signal reconstruction without requiring the sparsity level 𝐾. However, the BAOMP algorithm adopts two parameters (atom-adding constant 𝜇₁ and atom-deleting constant 𝜇₂) which are related with sparsity level and required to be tuned online.

In this paper, a novel method named as splitting matching pursuit (SMP) is proposed to reconstruct sparse signal in compressed sensing. The SP/CoSaMP algorithm assumes that the true indices in the support set of the original signal correspond to the large components of the correlation vector and choose 𝐾/2𝐾 largest components at each iteration. However, in practice some true indices may correspond to a set of small components. The proposed method selects 𝐹𝑙 (𝐹𝑙>2𝐾) largest components of the correlation vector, which expands the searching area and increases the probability of finding the true components at one iteration. The 𝐹𝑙 candidate components are divided into 𝐹 split sets with equal length 𝑙. The proposed method does not require the sparsity level 𝐾 to be known in prior. Different from BAOMP algorithm which adopts two tuning parameters related with sparsity level, the proposed algorithm uses two parameters 𝐹 and 𝑙 which are preset and determined by 𝑁 and 𝑀. The candidate components are divided into 𝐹 split sets, which could be processed simultaneously and suitable for parallel computation.

The proposed SMP method is then extended to more practical condition, for example, the direction of arrival (DOA) estimation problem in phased-array radar system using compressed sensing. Numerical simulations show that the proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods. The proposed method also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon [10], APES [11], and GLRT [12] based on hundreds of snapshots.

The rest of the sections are organized as follows. The proposed SMP method is introduced in Section 2. The simulation results are listed in Section 3, and the paper is summarized in Section 4.

2. The Splitting Matching Pursuit Method

The most difficult part of signal reconstruction is to identify the locations of 𝐾 largest components in the target signal. The SP/CoSaMP algorithm assumes that the true indices in the support set of the original signal correspond to the large components of the correlation vector and choose 𝐾/2𝐾 largest components at each iteration. However, in practice some true indices may correspond to a set of small components. The proposed SMP selects 𝐹𝑙 (𝐹𝑙>2𝐾) largest components of the correlation vector, which expands the searching area and increases the probability of finding the true components at one iteration. In the following, a short introduction of the proposed method is given in Section 2.1, and the detailed algorithm is introduced in Section 2.2. The convergency analysis, parameters setting, complexity analysis, and convergency speed analysis of the proposed method are provided in Sections 2.3 to 2.6, respectively.

2.1. Brief Introduction to Splitting Matching Pursuit Method

A schematic diagram of the proposed algorithm is depicted in Figure 1. At the beginning of each iteration, the proposed method selects 𝐹𝑙 (𝐹𝑙>2𝐾) largest components of the correlation vector, which are divided into 𝐹 split sets with equal length 𝑙. Each split set is merged with the estimated support set from the previous iteration, resulting in a sequence of merged split sets. The proposed algorithm then approximates the target signal on each merged split set to obtain a sequence of split estimates using least-square algorithm. A set of pruned split sets are then built by retaining only the 𝑙 largest magnitude entries in the split estimates. The obtained pruned split sets are then combined to a final merged set, which contains as much as possible true components. An interim estimate is then calculated based on the final merged set using least-square algorithm. An estimated support set is obtained by retaining 𝑙 indices corresponding to the largest magnitude entries in the interim estimate, based on which a final estimate is generated. The iterations repeat if the 𝑙₂ norm (magnitude) of the calculated residual is less than a threshold 𝑇.

Figure 1 

Description of reconstruction procedures of the SMP method.

2.2. Detailed Procedures of Splitting Matching Pursuit Method

The detailed procedure of the SMP method is listed as follows.

Algorithm 1 (the SMP method). Input. Sensing matrix Φ, measurement vector 𝑦, parameters 𝐹 and 𝑙, and threshold 𝑇 to halt the iterations.

Output. The estimated signal 𝑥_out.

Initialization(1)𝑎⁰=∅, where 𝑎⁰ indicates the initial estimated support set and ∅ denotes empty set. 𝑎 is the estimated support set with length 𝑙. (2)𝑟⁰=𝑦, where 𝑟⁰ denotes the initial residual. (3)𝑐⁰=Φ^∗𝑟⁰, where 𝑐⁰ denotes the initial correlation vector and Φ^∗ denotes the transpose of matrix Φ.

Iteration. At the 𝑘th iteration, go through the following steps.(1)Splitting. Locate the 𝐹𝑙 largest components of the correlation vector at the (𝑘−1)th iteration, and divide them into 𝐹 split sets as 𝐼^𝑘_𝑗=𝑐^𝑘−1_{(𝑗−1)𝑙+1∶𝑗𝑙},𝑗=1,…,𝐹, (6) where 𝐼^𝑘_𝑗 denotes the 𝑗th split set at the 𝑘th iteration, and 𝑐^𝑘−1 denotes the correlation vector at the (𝑘−1)th iteration. Furthermore, 𝑐^𝑘−1_{(𝑗−1)𝑙+1∶𝑗𝑙} denotes the ((𝑗−1)𝑙+1)th largest magnitude entry to the (𝑗𝑙)th largest magnitude entry of 𝑐^𝑘−1,𝑗=1,…,𝐹.(2)Support Merging. Each newly identified split set is united with the estimated support set from the previous iteration, resulting in a sequence of merged split sets as 𝐽^𝑘_𝑗=𝐼^𝑘_𝑗∪𝑎^𝑘−1,𝑗=1,…,𝐹, (7) where 𝐽^𝑘_𝑗 denotes the 𝑗th merged split set at the 𝑘th iteration, and 𝑎^𝑘−1 denotes the estimated support set at the (𝑘−1)th iteration.(3)Estimation. The proposed algorithm then solves a least square problem to approximate the nonzero entries of the target signal on each merged split set (𝐽^𝑘_𝑗,𝑗=1,…,𝐹) and sets other entries as zero, resulting in a sequence of split estimates as 𝑏^𝑘_{𝑗𝐽𝑘𝑗}=Φ_𝐽^𝑘_𝑗^†𝑏𝑦,_{𝑘𝑗(𝐽𝑘𝑗)𝑐}=0,𝑗=1,…,𝐹, (8) where 𝑏^𝑘_𝑗 denotes the 𝑗th split estimate of the target signal at the 𝑘th iteration. The vector 𝑏^𝑘_𝑗|_𝐽^𝑘_𝑗 is composed of the entries of   𝑏^𝑘_𝑗 indexed by 𝑖∈𝐽^𝑘_𝑗,   and 𝑏^𝑘_𝑗|_(𝐽^𝑘_𝑗)^𝑐 is composed of the entries of  𝑏^𝑘_𝑗 indexed by 𝑖∈(𝐽^𝑘_𝑗)^𝑐. † indicates pseudoinverse operation. The matrix Φ_𝐽^𝑘_𝑗 consists of the columns of Φ with indices 𝑖∈𝐽^𝑘_𝑗. (4)Pruning. Obtain a sequence of pruned split sets via retaining only the largest 𝑙 indices corresponding to the largest magnitude entries in 𝑏^𝑘_𝑗|_𝐽^𝑘_𝑗,𝑗=1,…,𝐹, for example, 𝐻^𝑘_𝑗=indicesof𝑙largestmagnitudeentriesin𝑏^𝑘_𝑗|_𝐽^𝑘_𝑗,𝑗=1,…,𝐹, (9) where 𝐻^𝑘_𝑗 denotes the 𝑗th pruned split set at the 𝑘th iteration. (5)Split Sets Merging. The pruned split sets are merged to form a final merged set, 𝐺^𝑘, as 𝐺^𝑘=union𝐻^𝑘₁,𝐻^𝑘₂,…,𝐻^𝑘_𝐹. (10)(6)Estimation. An interim estimate of the original signal is calculated based on the final merged set 𝐺^𝑘 using least-square algorithm as 𝑞^𝑘_𝐺𝑘=𝐺^𝑘†𝑞𝑦,_{𝑘(𝐺𝑘)𝑐}=0, (11) where 𝑞^𝑘 denotes the interim estimate of the original signal at the 𝑘th iteration. The vector 𝑞^𝑘|_𝐺^𝑘 is composed of the entries of 𝑞^𝑘 indexed by 𝑖∈𝐺^𝑘, and 𝑞^𝑘|_(𝐺^𝑘₎^𝑐 is composed of the entries of 𝑞^𝑘 indexed by 𝑖∈(𝐺^𝑘)^𝑐.(7)Pruning. Obtain the estimated support set by retaining 𝑙 indices corresponding to the largest magnitude entries in the vector 𝑞^𝑘|_𝐺^𝑘, as 𝑎^𝑘=indicesof𝑙largestmagnitudeentriesin𝑞^𝑘_𝐺𝑘, (12) where 𝑎^𝑘 denotes the estimated support set at the 𝑘th iteration. (8)Estimation. Obtain the final estimate at each iteration, based on 𝑎^𝑘 using least-square algorithm as 𝑥^𝑘_𝐹𝑎𝑘=Φ_𝑎𝑘^†𝑥𝑦,_{𝑘𝐹(𝑎𝑘)𝑐}=0, (13) where 𝑥^𝑘_𝐹 denotes the final estimate at the 𝑘th iteration. The vector 𝑥^𝑘_𝐹|_𝑎𝑘 is composed of the entries of 𝑥^𝑘_𝐹 indexed by 𝑖∈𝑎^𝑘, and 𝑥^𝑘_𝐹|_(𝑎𝑘)^𝑐 is composed of the entries of 𝑥^𝑘_𝐹 indexed by 𝑖∈(𝑎^𝑘)^𝑐. The matrix Φ_𝑎𝑘 consists of the columns of Φ with indices 𝑖∈𝑎^𝑘.(9) Residual calculation: 𝑟^𝑘=𝑦−Φ_𝑎𝑘𝑥^𝑘_𝐹𝑎𝑘, (14) where 𝑟^𝑘 denotes the residual at the 𝑘th iteration. (10) If ‖𝑟^𝑘‖₂>𝑇, perform the correlation calculation 𝑐^𝑘=Φ^∗𝑟^𝑘, and then go to step (1) of the (𝑘+1)th iteration; otherwise, set 𝑥_out=𝑥^𝑘_𝐹 and quit the iteration.

2.3. Convergency Analysis

Here, we will discuss the convergency of the proposed SMP method.

Theorem 2 (Theorem  2.1 in [8]). Let 𝑥∈ℜ^𝑁 be a 𝐾-sparse signal, and let its corresponding measurement be 𝑦=Φ𝑥+𝑒∈ℜ^𝑀. If the sampling matrix satisfies the restricted isometry property (RIP) with constant 𝛿_4𝐾<0.1, (15)then the signal approximation 𝑥^𝑘 is 𝐾-sparse and𝑥−𝑥^𝑘+1₂≤0.5𝑥−𝑥^𝑘₂+10𝑣. (16)
In particular, 𝑥−𝑥^𝑘₂≤2^−𝑘‖𝑥‖₂+20𝑣, (17) where 𝑣 denotes the unrecoverable energy in the signal.

Proposition 3. Let 𝑥∈ℜ^𝑁 be a 𝐾-sparse signal, and let its corresponding measurement be 𝑦=Φ𝑥+𝑒∈ℜ^𝑀. If the sampling matrix satisfies the RIP with constant 𝛿_𝐹𝑙<0.1, (18) then the proposed SMP algorithm is guaranteed to recover 𝑥 from 𝑦 via a finite number of iterations.

Proof. The proving process is very similar to that of Theorem 2 (Theorem  2.1 in [8]). The CoSaMP algorithm selects the 2𝐾 largest components of the correlation vector at each iteration, while the proposed SMP method selects the 𝐹𝑙 largest components of the correlation vector. We can obtain the similar results through replacing 2𝐾 with 𝐹𝑙 in the derivation process in [8].

2.4. Parameters Setting

Here, we will discuss how to set values of the number of split sets 𝐹 and the length of the split set 𝑙.

Proposition 4. Note that 𝐹𝑙<4𝐾 guarantees 𝛿_𝐹𝑙<0.1 if 𝛿_4𝐾<0.1.

Proof. Considering the monotonicity of 𝛿_𝐾: for any two integers 𝐾≤𝐾^′, 𝛿_𝐾≤𝛿^′_𝐾 according to Lemma  1 in [7]. So we have 𝛿_𝐹𝑙<𝛿_4𝐾 provided that 𝐹𝑙<4𝐾.

Proposition 5. Step  (7) in Algorithm 1 guarantees 𝐾≤𝑙.

Proof. As the estimated support set, 𝑎^𝑘 contains at least 𝐾 elements, resulting in 𝐾≤𝑙.
According to Propositions 3~5, in order to guarantee a perfect recovery of the 𝐾-sparse vector 𝑥 from 𝑦 via a finite number of iterations, we have 𝐹𝑙<4𝐾≤4𝑙, resulting in 𝐹<4. Since 𝐹 is set large enough to expand the searching area, 𝐹 is set as 3 in the proposed method. We then have 𝑙<(4/3)𝐾 according to 𝐹𝑙<4𝐾. Propositions 3~5 are based on Theorem 2, which has a rigid setting for 𝛿_4𝐾. We can relax the range of 𝑙 to [𝐾,2𝐾). As a result, we need not know the exact value of 𝐾 and can select an integer randomly from [𝐾,2𝐾) based on an estimated value of 𝐾.

2.5. Complexity Analysis

The process of complexity analysis is similar to that of [7]. In each iteration, the correlation maximization procedure requires 𝑀𝑁 computations in general, while the cost of computing the projections is of the order of 𝑂(𝐹(2𝑠)²𝑀+(2𝐹𝑠)²𝑀), which could be approximated as 𝑂(𝑠²𝑀) when 𝐹 is chosen as a small value (e.g. 3 in the simulation setup). As a result, the total cost of computing is of the order of 𝑂(𝑀𝑁+𝑠²𝑀), which is comparable to the SP/CoSaMP algorithm.

2.6. Convergence Speed Analysis

Here, we will discuss the convergence speed of the proposed SMP method.

Theorem 6 (Theorem  8 in [7]). The number of iterations (𝑛^𝑆𝑃_𝑖𝑡) of the SP algorithm is upper bounded by 𝑛^𝑆𝑃_𝑖𝑡≤1.5𝐾−log𝑐_𝐾, (19) where 𝑐_𝐾=2𝛿_3𝐾(1+𝛿_3𝐾)/(1−𝛿_3𝐾)³.

Proposition 7. The number of iterations (𝑛^𝑆𝑀𝑃_𝑖𝑡) of the SMP algorithm is upper bounded by 𝑛^𝑆𝑀𝑃_𝑖𝑡<𝑛^SP_𝑖𝑡𝐹. (20)

Proof. The SP algorithm selects the 𝐾 largest components of the correlation vector at each iteration, while the proposed SMP method selects the 𝐹𝑙 (𝐹𝑙≥𝐹𝐾 according to Proposition 5) largest components of the correlation vector. The searching area is thus expanded to at least 𝐹 times of that of SP algorithm. The probability of finding the true components at one iteration is increased to at least 𝐹 times of that of SP algorithm. As a result, the number of iterations is decreased to 1/𝐹 of that of the SP algorithm by average according to Theorem 6.

3. Simulation Results and Analysis

In this section, a simple example is firstly carried out to verify the performance of the proposed SMP method in reconstructing zero-one binary and Gaussian signals. The proposed method is then extended to cope with the DOA estimation problem in phased-array radar system.

3.1. A Simple Example

The simulations are carried out to compare the accuracy of different reconstruction algorithms empirically. The proposed SMP algorithm is compared with some popular greedy algorithms (including OMP, ROMP, StOMP, SP, CoSaMP, and BAOMP algorithms) and the convex optimization algorithm (BP method).

In the simulation setup, a signal sparsity level 𝐾 is chosen such that 𝐾≤𝑀/2 given the length of the original signal 𝑁 (𝑁=256) and the length of the measurement vector 𝑀 (𝑀=128). An 𝑀×𝑁 sampling matrix Φ is randomly generated from the standard i.i.d. Gaussian ensemble. A support set 𝑆 of size 𝐾 is selected uniformly at random, and the original sparse signal vector 𝑥 is chosen as either Gaussian signal or zero-one signal [7]. The estimate of the original signal, 𝑥, is computed based on the measurement vector 𝑦 generated through 𝑦=Φ𝑥. In the experiments, OMP uses 𝐾 iterations, StOMP and BP methods use the default settings (OMP, StOMP, and BP tools use SparseLab [13]), and ROMP and SP methods use the parameters given in [6, 7], respectively. The proposed SMP method use 𝐹=3, 𝑙=84, 𝑛_max=𝑀, 𝜀=10⁻⁵, as the input parameters.

The signal sparsity level 𝐾 is varied from 0 to 𝑀/2. For each fixed 𝐾, five hundred Monte Carlo simulations are carried out for each algorithm. The reconstruction is considered to be exact when the 𝑙₂ norm of the difference between the original signal 𝑥 and the reconstructed one 𝑥 is smaller than 10⁻⁵, that is, ‖𝑥−𝑥‖₂<10⁻⁵. The frequency of exact reconstruction (𝜉) is used to evaluate the reconstruction performance of the different methods, which is defined as 𝛼𝜉=𝑁_MC, (21) where 𝛼 denotes the number of exact reconstructions for each algorithm given a fixed 𝐾, and 𝑁_MC denotes the number of Monte Carlo simulations. The frequency of exact reconstruction is also adopted by the SP method to evaluate the reconstruction performance [7].

Figures 2 and 3 show the reconstruction results for binary zero-one and Gaussian sparse signals, respectively. We only present the results of the SP algorithm since the SP and CoSaMP algorithms are almost the same with different deviation process, and both obtain the same simulation results. As can be seen in Figure 2, for binary zero-one sparse signal which is a difficult case for OMP-type methods, the performance of the proposed SMP method is much better than all other OMP-type methods and comparable to the BP minimization method. Of particular interest is the sparsity level at which the recovery rate drops below 100%, that is, the critical sparsity defined in [7]. It could be seen from Figure 2 that the proposed method is with the largest critical sparsity (39), which exceeds that of the BP method (36). For the Gaussian sparse signal, as shown in Figure 3, the proposed method also gives the comparable performance to the BP method.

Figure 2 

Zero-one signal.

Figure 3 

Gaussian signal.

3.2. DOA Estimation Based on Splitting Matching Pursuit Method

In this section, the proposed SMP method is extended to solve the DOA estimation problem in phased-array radar system. The signal model for DOA estimation in phased-array radar system is represented in a standard compressed sensing form in Section 3.2.1, where the sparse radar scene is abstracted as a sparse signal. And the simulation results are shown in Section 3.2.2, which shows that the proposed method succeeds in identifying multiple targets in a sparse radar scene.

3.2.1. Signal Model for DOA Estimation and Sparse Representation

Assume that a phased-array radar system consists of half wavelength spaced uniform linear arrays (ULAs). Targets may appear at directions represented by DOA angles. The task of signal processing is to estimate the directions to the targets and the corresponding complex amplitudes (DOA estimation, see [14]). We assume that the other parameters like range and Doppler frequency have been isolated before by appropriate processing.

The ULA of the phased-array radar system consists of 𝑀 antennas, which are used to emit the transmitted signal 𝑠(𝑡). The 𝑀×1 received complex vector of array observations is defined as 𝑓(𝑡)=[𝑓₁(𝑡),…,𝑓_𝑀(𝑡)]^𝑇. Assuming a hypothetical target located at a DOA angle of 𝜃 in the far field, the received complex vector of array observations can be written as 𝑓(𝑡)=𝛽(𝜃)𝑠(𝑡)𝑎(𝜃)+𝑛(𝑡), (22) where 𝛽(𝜃) is the reflection coefficient of the hypothetical target, and 𝑛(𝑡) is an 𝑀×1 complex Gaussian noise vector. 𝑎(𝜃) is the 𝑀×1 steering vector, which is defined as 𝑎(𝜃)=1𝑒^{𝑗(2𝜋𝑑sin𝜃/𝜆)},…,𝑒^{𝑗(𝑀−1)(2𝜋𝑑sin𝜃/𝜆)𝑇}, (23) where 𝑑 is the distance between the elements of the arrays and 𝜆 denotes wavelength.

Assuming 𝐷 targets are observed with reflection coefficients {𝛽_𝑖}^𝐷_𝑖=1 and DOA angles {𝜃_𝑖}^𝐷_𝑖=1, the 𝑀×1 received complex vector of array observations can be written as 𝛼(𝑡)=^𝐷_𝑖=1𝛽𝜃_𝑖𝑠(𝑡)𝑎𝜃_𝑖+𝛾(𝑡), (24) where 𝛾(𝑡) is an 𝑀×1 complex Gaussian noise vector. Equation (24) could be rewritten as 𝛼(𝑡)=𝐴(𝜃)𝑆(𝑡,𝜃)𝛾(𝑡), (25) where 𝐴(𝜃)=[𝑎(𝜃₁)𝑎(𝜃₂)⋯𝑎(𝜃_𝐷)] is an 𝑀×𝐷 steering matrix, and 𝑆(𝑡,𝜃)=𝑠(𝑡)[𝛽(𝜃₁)𝛽(𝜃₂)⋯𝛽(𝜃_𝐷)]^𝑇 denotes a 𝐷×1 reflection vector.

Since the radar scene is generally in practice sparse, compressed sensing is a valid candidate for estimating the DOA angles for multiple targets. To do so, the DOA angle plane is divided into 𝑁 fine grids, each cell generally with the same size Δ𝜃. The 𝑖th grid represents the DOA angle of 𝜃₀+(𝑖−1)Δ𝜃, where 𝜃₀ is the initial angle of the DOA plane. The steering matrix and reflection vector in (25) are extended to obtain the 𝑀×𝑁 extended steering matrixΦ and the 𝑁×1 extended reflection vector 𝑥, which are defined as Φ=[𝑎(𝜃₀)𝑎(𝜃₀+Δ𝜃)⋯𝑎(𝜃₀+(𝑖−1)Δ𝜃)𝑎(𝜃₀+(𝑁−1)Δ𝜃)] and 𝑥=𝑠(𝑡)[𝛽(𝜃₀)𝛽(𝜃₀+Δ𝜃)⋯𝛽(𝜃₀+(𝑖−1)Δ𝜃)𝛽(𝜃₀+(𝑁−1)Δ𝜃)]^𝑇. Since small number of grids are occupied by the targets, 𝑥 is a sparse vector with the 𝑖th element defined as 𝑥(𝑖)=𝑠(𝑡)𝛽(𝜃₀+(𝑖−1)Δ𝜃) if the 𝑖th grid is occupied by the target; otherwise, 𝑥(𝑖)=0. As a result, the 𝑀×1 received complex vector of array observations 𝑦 could be written as 𝑦=Φ𝑥+𝑒, (26) where 𝑒 is an 𝑀×1 complex Gaussian noise vector. Though in (26) the radar vectors and matrices are complex valued in contrary to the original compressed sensing environment, it is easy to transfer it to real variables according to [15, 16].

Discussion. In [17, 18], it is assumed that the discretized step is small enough so that each target falls on some specific grid point. However, no matter how finely the parameter space is gridded, the sources may not lie in the center of the grid cells, and consequently there is mismatch between the assumed and the actual bases for sparsity. The sensitivity of compressed sensing to mismatch between the assumed and the actual sparsity bases is studied in [19]. The effect of basis mismatch is analyzed on the best 𝑘-term approximation error, and some achievable bounds for the 𝑙₁ error of the best 𝑘-term approximation are provided. The readers can refer to [19] for a detailed analysis on the influence of the griding operations on the estimation performance.

3.2.2. Simulation Results

In this section, an example about DOA estimation is provided based on the phased-array radar system, which consists of half wavelength spaced uniform linear arrays (ULAs). The number of transmit/receive antennas is 20. The antennas transmit independent orthogonal quadrature phase shift keyed (QPSK) waveforms and the carrier frequency is 8.62 GHz. The SNR of the measurement noise is set to a fixed value (20 dB). The range of the DOA plane is [0^∘,90^∘], which is divided into 30 cells with the initial angle (𝜃₀) and angle interval (Δ𝜃) equaling 0^∘ and 3^∘, respectively. A maximum of 𝐿=512 snapshots are considered at the receive node. Targets may appear at directions represented by DOA angles.

The proposed SMP method is compared to several popular greedy algorithms, for example, the OMP, ROMP, StOMP, SP, CoSaMP and BAOMP algorithms, and the convex optimization algorithm, BPDN method. The proposed method is further compared to three commonly used methods in DOA estimation, for example, the Capon, APES, and GLRT methods. The compressed sensing based methods (the SMP method, the greedy algorithms, and the BPDN method) use one snapshot only, and the Capon, APES, and GLRT methods use 512 snapshots each. Five hundred Monte Carlo simulations are carried out, and in each trial four targets locate randomly within the DOA range of [0^∘,90^∘], and the corresponding reflection coefficients are set as {𝛽_𝑘=1,𝑘=1,…,4}.

The average reconstruction error is adopted to evaluate the reconstruction performance of the methods, which is defined as 𝜒_average=^𝑁MC_𝑖=1𝜒^𝑖𝑁_MC, (27) where 𝜒^𝑖 denotes the reconstruction error at the 𝑖th Monte Carlo simulation, which is defined as 𝜒^𝑖=𝑥^𝑖_estimate−𝑥^𝑖₂𝑥^𝑖₂, (28) where 𝑥^𝑖 and 𝑥^𝑖_estimate represent the true and estimated signal representing the sparse radar scene at the 𝑖th Monte Carlo simulation, respectively. The average root mean square error (RMSE) is also adopted to evaluate the DOA estimation performance of the methods, which is defined in [20].

The results in Table 1 show that the proposed SMP method is with the smallest average reconstruction error and average RMSE. The proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods. It also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon, APES, and GLRT based on 512 snapshots.

Table 1 

Performance comparison.

4. Conclusion

We have presented a novel SMP method for sparse signal reconstruction in compressed sensing. The proposed method expands the searching area and increases the probability of finding the true components at one iteration. It also does not require the sparsity level 𝐾 to be known in prior. The proposed method is then extended to more practical condition, for example, the direction of arrival (DOA) estimation problem in phased-array radar system using compressed sensing. Numerical simulations show that the proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods. The proposed method also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon, APES, and GLRT based on hundreds of snapshots.

Acknowledgments

This work was supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (61221063), the Natural Science Foundations of China (no. 61104051, no. 61074176, and no. 61004087), the Fundamental Research Funds for the Central Universities, and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.