# Splitting Matching Pursuit Method for Reconstructing Sparse Signal in Compressed Sensing

Authors: Liu Jing; Han ChongZhao; Yao XiangHua; Lian Feng
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

In this paper, a novel method named as splitting matching pursuit (SMP) is proposed to reconstructK-sparse
signal in compressed sensing. The proposed method selects Fl  (Fl>2K) largest components of the correlation
vector c, which are divided into F split sets with equal length l. 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 K 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.


---

## Body

## 1. Introduction

The standard noiseless model in compressed sensing is(1)y=Φx,
where x∈ℜN is a K-sparse signal (K≪N), y∈ℜM is a measurement of x, and Φ is an M×N sensing matrix. The compressed sensing recovery problem is defined as follows: given y and Φ, find a signal x within the class of interest satisfies (1) exactly. The compressed sensing recovery process consists of a search for the sparsest signal x that yields the measurement y. By defining the l0 “norm” of a vector ∥w∥0 as the number of nonzero entries in w, the simplest way to pose a recovery algorithm is using the optimization(2)x=argminy=Φw∥w∥0.Solutions to (2) therefore lead to algorithms for recovering K-sparse signals from M linear measurements. In general, the minimization of (2) is NP-hard. An alternative to the l0 “norm” used in (2) is to use the l1 “norm”, defined as ∥w∥1=∑n=1N|w(n)|. The resulting adaptation of (2), known as basis pursuit (BP) [1], is formally defined as(3)x=argminy=Φw∥w∥1.Since thel1 “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 y=Φx+e, where e denotes an M×1 measurement noise vector. We simply change the constraint on the solution to
(4)x=argmin∥y-Φw∥2≤ϵ∥w∥1,
where ϵ>∥e∥2 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
(5)x=argmin∥w∥1+λ∥y-Φw∥2
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 K. However, the BAOMP algorithm adopts two parameters (atom-adding constant μ1 and atom-deleting constant μ2) 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 chooseK/2K largest components at each iteration. However, in practice some true indices may correspond to a set of small components. The proposed method selects Fl (Fl>2K) largest components of the correlation vector, which expands the searching area and increases the probability of finding the true components at one iteration. The Fl candidate components are divided into F split sets with equal length l. The proposed method does not require the sparsity level K to be known in prior. Different from BAOMP algorithm which adopts two tuning parameters related with sparsity level, the proposed algorithm uses two parameters F and l which are preset and determined by N and M. The candidate components are divided into F 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 Section2. 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 ofK 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 K/2K largest components at each iteration. However, in practice some true indices may correspond to a set of small components. The proposed SMP selects Fl (Fl>2K) 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 Figure1. At the beginning of each iteration, the proposed method selects Fl (Fl>2K) largest components of the correlation vector, which are divided into F split sets with equal length l. 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 l 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 l indices corresponding to the largest magnitude entries in the interim estimate, based on which a final estimate is generated. The iterations repeat if the l2 norm (magnitude) of the calculated residual is less than a threshold T.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 y, parameters F and l, and threshold T to halt the iterations.Output. The estimated signal xout.Initialization
(1)
a
0
=
∅, where a0 indicates the initial estimated support set and ∅ denotes empty set. a is the estimated support set with length l.
(2)
r
0
=
y, where r0 denotes the initial residual.
(3)
c
0
=
Φ
*
r
0, where c0 denotes the initial correlation vector and Φ* denotes the transpose of matrix Φ.Iteration. At the kth iteration, go through the following steps.(1)
Splitting. Locate the Fl largest components of the correlation vector at the (k-1)th iteration, and divide them into F split sets as
(6)Ijk=c(j-1)l+1:jlk-1,j=1,…,F,
where Ijk denotes the jth split set at the kth iteration, and ck-1 denotes the correlation vector at the (k-1)th iteration. Furthermore, c(j-1)l+1:jlk-1 denotes the ((j-1)l+1)th largest magnitude entry to the (jl)th largest magnitude entry of ck-1,j=1,…,F.(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
(7)Jjk=Ijk∪ak-1,j=1,…,F,
where Jjk denotes the jth merged split set at the kth iteration, and ak-1 denotes the estimated support set at the (k-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 (Jjk,j=1,…,F) and sets other entries as zero, resulting in a sequence of split estimates as
(8)bjk|Jjk=(ΦJjk)†y,bjk|(Jjk)c=0,j=1,…,F,
where bjk denotes the jth split estimate of the target signal at the kth iteration. The vector bjk|Jjk is composed of the entries of   bjk indexed by i∈Jjk,   and bjk|(Jjk)c is composed of the entries of  bjk indexed by i∈(Jjk)c. † indicates pseudoinverse operation. The matrix ΦJjk consists of the columns of Φ with indices i∈Jjk.(4)
Pruning. Obtain a sequence of pruned split sets via retaining only the largest l indices corresponding to the largest magnitude entries in bjk|Jjk,j=1,…,F, for example,
(9)Hjk={indicesofllargestmagnitudeentriesinbjk|Jjk},j=1,…,F,
where Hjk denotes the jth pruned split set at the kth iteration.(5)
Split Sets Merging. The pruned split sets are merged to form a final merged set, Gk, as
(10)Gk=union{H1k,H2k,…,HFk}.(6)
Estimation. An interim estimate of the original signal is calculated based on the final merged set Gk using least-square algorithm as
(11)qk|Gk=(Gk)†y,qk|(Gk)c=0,
where qk denotes the interim estimate of the original signal at the kth iteration. The vector qk|Gk is composed of the entries of qk indexed by i∈Gk, and qk|(Gk)c is composed of the entries of qk indexed by i∈(Gk)c.(7)
Pruning. Obtain the estimated support set by retaining l indices corresponding to the largest magnitude entries in the vector qk|Gk, as
(12)ak={indicesofllargestmagnitudeentriesinqk|Gk},
where ak denotes the estimated support set at the kth iteration.(8)
Estimation. Obtain the final estimate at each iteration, based on ak using least-square algorithm as
(13)xFk|ak=(Φak)†y,xFk|(ak)c=0,
where xFk denotes the final estimate at the kth iteration. The vector xFk|ak is composed of the entries of xFk indexed by i∈ak, and xFk|(ak)c is composed of the entries of xFk indexed by i∈(ak)c. The matrix Φak consists of the columns of Φ with indices i∈ak.(9)
Residual calculation:(14)rk=y-(Φak)xFk|ak,
where rk denotes the residual at the kth iteration.(10)
If∥rk∥2>T, perform the correlation calculation ck=Φ*rk, and then go to step (1) of the (k+1)th iteration; otherwise, set xout=xFk 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]).
Letx∈ℜN be a K-sparse signal, and let its corresponding measurement be y=Φx+e∈ℜM. If the sampling matrix satisfies the restricted isometry property (RIP) with constant
(15)δ4K<0.1,then the signal approximation x^k is K-sparse and(16)∥x-x^k+1∥2≤0.5∥x-x^k∥2+10v.
In particular,(17)∥x-x^k∥2≤2-k∥x∥2+20v,
where v denotes the unrecoverable energy in the signal.Proposition 3.
Letx∈ℜN be a K-sparse signal, and let its corresponding measurement be y=Φx+e∈ℜM. If the sampling matrix satisfies the RIP with constant
(18)δFl<0.1,
then the proposed SMP algorithm is guaranteed to recover x from y via a finite number of iterations.Proof.
The proving process is very similar to that of Theorem2 (Theorem  2.1 in [8]). The CoSaMP algorithm selects the 2K largest components of the correlation vector at each iteration, while the proposed SMP method selects the Fl largest components of the correlation vector. We can obtain the similar results through replacing 2K with Fl in the derivation process in [8].


### 2.4. Parameters Setting

Here, we will discuss how to set values of the number of split setsF and the length of the split set l.Proposition 4.
Note thatFl<4K guarantees δFl<0.1 if δ4K<0.1.Proof.
Considering the monotonicity ofδK: for any two integers K≤K′, δK≤δK′ according to Lemma  1 in [7]. So we have δFl<δ4K provided that Fl<4K.Proposition 5.
Step(7) in Algorithm 1 guarantees K≤l.Proof.
As the estimated support set,ak contains at least K elements, resulting in K≤l.
According to Propositions3~5, in order to guarantee a perfect recovery of the K-sparse vector x from y via a finite number of iterations, we have Fl<4K≤4l, resulting in F<4. Since F is set large enough to expand the searching area, F is set as 3 in the proposed method. We then have l<(4/3)K according to Fl<4K. Propositions 3~5 are based on Theorem 2, which has a rigid setting for δ4K. We can relax the range of l to [K,2K). As a result, we need not know the exact value of K and can select an integer randomly from [K,2K) based on an estimated value of K.


### 2.5. Complexity Analysis

The process of complexity analysis is similar to that of [7]. In each iteration, the correlation maximization procedure requires MN computations in general, while the cost of computing the projections is of the order of O(F(2s)2M+(2Fs)2M), which could be approximated as O(s2M) when F 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 O(MN+s2M), 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(nitSP) of the SP algorithm is upper bounded by
(19)nitSP≤1.5K-log(cK),
where cK=2δ3K(1+δ3K)/(1-δ3K)3.Proposition 7.
The number of iterations(nitSMP) of the SMP algorithm is upper bounded by
(20)nitSMP<nitSPF.Proof.
The SP algorithm selects theK largest components of the correlation vector at each iteration, while the proposed SMP method selects the Fl (Fl≥FK according to Proposition 5) largest components of the correlation vector. The searching area is thus expanded to at least F times of that of SP algorithm. The probability of finding the true components at one iteration is increased to at least F times of that of SP algorithm. As a result, the number of iterations is decreased to 1/F of that of the SP algorithm by average according to Theorem 6.


## 2.1. Brief Introduction to Splitting Matching Pursuit Method

A schematic diagram of the proposed algorithm is depicted in Figure1. At the beginning of each iteration, the proposed method selects Fl (Fl>2K) largest components of the correlation vector, which are divided into F split sets with equal length l. 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 l 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 l indices corresponding to the largest magnitude entries in the interim estimate, based on which a final estimate is generated. The iterations repeat if the l2 norm (magnitude) of the calculated residual is less than a threshold T.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 y, parameters F and l, and threshold T to halt the iterations.Output. The estimated signal xout.Initialization
(1)
a
0
=
∅, where a0 indicates the initial estimated support set and ∅ denotes empty set. a is the estimated support set with length l.
(2)
r
0
=
y, where r0 denotes the initial residual.
(3)
c
0
=
Φ
*
r
0, where c0 denotes the initial correlation vector and Φ* denotes the transpose of matrix Φ.Iteration. At the kth iteration, go through the following steps.(1)
Splitting. Locate the Fl largest components of the correlation vector at the (k-1)th iteration, and divide them into F split sets as
(6)Ijk=c(j-1)l+1:jlk-1,j=1,…,F,
where Ijk denotes the jth split set at the kth iteration, and ck-1 denotes the correlation vector at the (k-1)th iteration. Furthermore, c(j-1)l+1:jlk-1 denotes the ((j-1)l+1)th largest magnitude entry to the (jl)th largest magnitude entry of ck-1,j=1,…,F.(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
(7)Jjk=Ijk∪ak-1,j=1,…,F,
where Jjk denotes the jth merged split set at the kth iteration, and ak-1 denotes the estimated support set at the (k-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 (Jjk,j=1,…,F) and sets other entries as zero, resulting in a sequence of split estimates as
(8)bjk|Jjk=(ΦJjk)†y,bjk|(Jjk)c=0,j=1,…,F,
where bjk denotes the jth split estimate of the target signal at the kth iteration. The vector bjk|Jjk is composed of the entries of   bjk indexed by i∈Jjk,   and bjk|(Jjk)c is composed of the entries of  bjk indexed by i∈(Jjk)c. † indicates pseudoinverse operation. The matrix ΦJjk consists of the columns of Φ with indices i∈Jjk.(4)
Pruning. Obtain a sequence of pruned split sets via retaining only the largest l indices corresponding to the largest magnitude entries in bjk|Jjk,j=1,…,F, for example,
(9)Hjk={indicesofllargestmagnitudeentriesinbjk|Jjk},j=1,…,F,
where Hjk denotes the jth pruned split set at the kth iteration.(5)
Split Sets Merging. The pruned split sets are merged to form a final merged set, Gk, as
(10)Gk=union{H1k,H2k,…,HFk}.(6)
Estimation. An interim estimate of the original signal is calculated based on the final merged set Gk using least-square algorithm as
(11)qk|Gk=(Gk)†y,qk|(Gk)c=0,
where qk denotes the interim estimate of the original signal at the kth iteration. The vector qk|Gk is composed of the entries of qk indexed by i∈Gk, and qk|(Gk)c is composed of the entries of qk indexed by i∈(Gk)c.(7)
Pruning. Obtain the estimated support set by retaining l indices corresponding to the largest magnitude entries in the vector qk|Gk, as
(12)ak={indicesofllargestmagnitudeentriesinqk|Gk},
where ak denotes the estimated support set at the kth iteration.(8)
Estimation. Obtain the final estimate at each iteration, based on ak using least-square algorithm as
(13)xFk|ak=(Φak)†y,xFk|(ak)c=0,
where xFk denotes the final estimate at the kth iteration. The vector xFk|ak is composed of the entries of xFk indexed by i∈ak, and xFk|(ak)c is composed of the entries of xFk indexed by i∈(ak)c. The matrix Φak consists of the columns of Φ with indices i∈ak.(9)
Residual calculation:(14)rk=y-(Φak)xFk|ak,
where rk denotes the residual at the kth iteration.(10)
If∥rk∥2>T, perform the correlation calculation ck=Φ*rk, and then go to step (1) of the (k+1)th iteration; otherwise, set xout=xFk 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]).
Letx∈ℜN be a K-sparse signal, and let its corresponding measurement be y=Φx+e∈ℜM. If the sampling matrix satisfies the restricted isometry property (RIP) with constant
(15)δ4K<0.1,then the signal approximation x^k is K-sparse and(16)∥x-x^k+1∥2≤0.5∥x-x^k∥2+10v.
In particular,(17)∥x-x^k∥2≤2-k∥x∥2+20v,
where v denotes the unrecoverable energy in the signal.Proposition 3.
Letx∈ℜN be a K-sparse signal, and let its corresponding measurement be y=Φx+e∈ℜM. If the sampling matrix satisfies the RIP with constant
(18)δFl<0.1,
then the proposed SMP algorithm is guaranteed to recover x from y via a finite number of iterations.Proof.
The proving process is very similar to that of Theorem2 (Theorem  2.1 in [8]). The CoSaMP algorithm selects the 2K largest components of the correlation vector at each iteration, while the proposed SMP method selects the Fl largest components of the correlation vector. We can obtain the similar results through replacing 2K with Fl in the derivation process in [8].


## 2.4. Parameters Setting

Here, we will discuss how to set values of the number of split setsF and the length of the split set l.Proposition 4.
Note thatFl<4K guarantees δFl<0.1 if δ4K<0.1.Proof.
Considering the monotonicity ofδK: for any two integers K≤K′, δK≤δK′ according to Lemma  1 in [7]. So we have δFl<δ4K provided that Fl<4K.Proposition 5.
Step(7) in Algorithm 1 guarantees K≤l.Proof.
As the estimated support set,ak contains at least K elements, resulting in K≤l.
According to Propositions3~5, in order to guarantee a perfect recovery of the K-sparse vector x from y via a finite number of iterations, we have Fl<4K≤4l, resulting in F<4. Since F is set large enough to expand the searching area, F is set as 3 in the proposed method. We then have l<(4/3)K according to Fl<4K. Propositions 3~5 are based on Theorem 2, which has a rigid setting for δ4K. We can relax the range of l to [K,2K). As a result, we need not know the exact value of K and can select an integer randomly from [K,2K) based on an estimated value of K.


## 2.5. Complexity Analysis

The process of complexity analysis is similar to that of [7]. In each iteration, the correlation maximization procedure requires MN computations in general, while the cost of computing the projections is of the order of O(F(2s)2M+(2Fs)2M), which could be approximated as O(s2M) when F 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 O(MN+s2M), 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(nitSP) of the SP algorithm is upper bounded by
(19)nitSP≤1.5K-log(cK),
where cK=2δ3K(1+δ3K)/(1-δ3K)3.Proposition 7.
The number of iterations(nitSMP) of the SMP algorithm is upper bounded by
(20)nitSMP<nitSPF.Proof.
The SP algorithm selects theK largest components of the correlation vector at each iteration, while the proposed SMP method selects the Fl (Fl≥FK according to Proposition 5) largest components of the correlation vector. The searching area is thus expanded to at least F times of that of SP algorithm. The probability of finding the true components at one iteration is increased to at least F times of that of SP algorithm. As a result, the number of iterations is decreased to 1/F 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 levelK is chosen such that K≤M/2 given the length of the original signal N (N=256) and the length of the measurement vector M (M=128). An M×N sampling matrix Φ is randomly generated from the standard i.i.d. Gaussian ensemble. A support set S of size K is selected uniformly at random, and the original sparse signal vector x is chosen as either Gaussian signal or zero-one signal [7]. The estimate of the original signal, x^, is computed based on the measurement vector y generated through y=Φx. In the experiments, OMP uses K 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 F=3, l=84, nmax=M, ε=10-5, as the input parameters.The signal sparsity levelK is varied from 0 to M/2. For each fixed K, five hundred Monte Carlo simulations are carried out for each algorithm. The reconstruction is considered to be exact when the l2 norm of the difference between the original signal x and the reconstructed one x^ is smaller than 10-5, that is, ∥x-x^∥2<10-5. The frequency of exact reconstruction (ξ) is used to evaluate the reconstruction performance of the different methods, which is defined as
(21)ξ=αNMC,
where α denotes the number of exact reconstructions for each algorithm given a fixed K, and NMC 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].Figures2 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 Section3.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 ofM antennas, which are used to emit the transmitted signal s(t). The M×1 received complex vector of array observations is defined as f(t)=[f1(t),…,fM(t)]T. 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)f(t)=β(θ)s(t)a(θ)+n(t),
where β(θ) is the reflection coefficient of the hypothetical target, and n(t) is an M×1 complex Gaussian noise vector. a(θ) is the M×1 steering vector, which is defined as
(23)a(θ)=[1ej(2πdsinθ/λ),…,ej(M-1)(2πdsinθ/λ)]T,
where d is the distance between the elements of the arrays and λ denotes wavelength.AssumingD targets are observed with reflection coefficients {βi}i=1D and DOA angles {θi}i=1D, the M×1 received complex vector of array observations can be written as
(24)α(t)=∑i=1Dβ(θi)s(t)a(θi)+γ(t),
where γ(t) is an M×1 complex Gaussian noise vector. Equation (24) could be rewritten as
(25)α(t)=A(θ)S(t,θ)γ(t),
where A(θ)=[a(θ1)a(θ2)⋯a(θD)] is an M×D steering matrix, and S(t,θ)=s(t)[β(θ1)β(θ2)⋯β(θD)]T denotes a D×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 intoN fine grids, each cell generally with the same size Δθ. The ith grid represents the DOA angle of θ0+(i-1)Δθ, where θ0 is the initial angle of the DOA plane. The steering matrix and reflection vector in (25) are extended to obtain the M×N extended steering matrixΦ and the N×1 extended reflection vector x, which are defined as Φ=[a(θ0)a(θ0+Δθ)⋯a(θ0+(i-1)Δθ)a(θ0+(N-1)Δθ)] and x=s(t)[β(θ0)β(θ0+Δθ)⋯β(θ0+(i-1)Δθ)β(θ0+(N-1)Δθ)]T. Since small number of grids are occupied by the targets, x is a sparse vector with the ith element defined as x(i)=s(t)β(θ0+(i-1)Δθ) if the ith grid is occupied by the target; otherwise, x(i)=0. As a result, the M×1 received complex vector of array observations y could be written as
(26)y=Φx+e,
where e is an M×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 k-term approximation error, and some achievable bounds for the l1 error of the best k-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 is20. 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 (θ0) and angle interval (Δθ) equaling 0° and 3°, respectively. A maximum of L=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 use512 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 {βk=1,k=1,…,4}.The average reconstruction error is adopted to evaluate the reconstruction performance of the methods, which is defined as(27)χaverage=∑i=1NMCχiNMC,
where χi denotes the reconstruction error at the ith Monte Carlo simulation, which is defined as
(28)χi=∥xestimatei-xi∥2∥xi∥2,
where xi and xestimatei represent the true and estimated signal representing the sparse radar scene at the ith 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 Table1 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.
Average reconstruction error
Average RMSE of DOA angles (degree)
SMP
0.04
0.05
OMP
0.2
0.19
ROMP
0.24
0.23
StOMP
0.23
0.25
SP/CoSaMP
0.15
0.13
BAOMP
0.17
0.16
BPDN
0.09
0.07
Capon
0.38
1.1
APES
0.35
0.54
GLRT
0.27
0.17


## 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 levelK is chosen such that K≤M/2 given the length of the original signal N (N=256) and the length of the measurement vector M (M=128). An M×N sampling matrix Φ is randomly generated from the standard i.i.d. Gaussian ensemble. A support set S of size K is selected uniformly at random, and the original sparse signal vector x is chosen as either Gaussian signal or zero-one signal [7]. The estimate of the original signal, x^, is computed based on the measurement vector y generated through y=Φx. In the experiments, OMP uses K 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 F=3, l=84, nmax=M, ε=10-5, as the input parameters.The signal sparsity levelK is varied from 0 to M/2. For each fixed K, five hundred Monte Carlo simulations are carried out for each algorithm. The reconstruction is considered to be exact when the l2 norm of the difference between the original signal x and the reconstructed one x^ is smaller than 10-5, that is, ∥x-x^∥2<10-5. The frequency of exact reconstruction (ξ) is used to evaluate the reconstruction performance of the different methods, which is defined as
(21)ξ=αNMC,
where α denotes the number of exact reconstructions for each algorithm given a fixed K, and NMC 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].Figures2 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 Section3.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 ofM antennas, which are used to emit the transmitted signal s(t). The M×1 received complex vector of array observations is defined as f(t)=[f1(t),…,fM(t)]T. 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)f(t)=β(θ)s(t)a(θ)+n(t),
where β(θ) is the reflection coefficient of the hypothetical target, and n(t) is an M×1 complex Gaussian noise vector. a(θ) is the M×1 steering vector, which is defined as
(23)a(θ)=[1ej(2πdsinθ/λ),…,ej(M-1)(2πdsinθ/λ)]T,
where d is the distance between the elements of the arrays and λ denotes wavelength.AssumingD targets are observed with reflection coefficients {βi}i=1D and DOA angles {θi}i=1D, the M×1 received complex vector of array observations can be written as
(24)α(t)=∑i=1Dβ(θi)s(t)a(θi)+γ(t),
where γ(t) is an M×1 complex Gaussian noise vector. Equation (24) could be rewritten as
(25)α(t)=A(θ)S(t,θ)γ(t),
where A(θ)=[a(θ1)a(θ2)⋯a(θD)] is an M×D steering matrix, and S(t,θ)=s(t)[β(θ1)β(θ2)⋯β(θD)]T denotes a D×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 intoN fine grids, each cell generally with the same size Δθ. The ith grid represents the DOA angle of θ0+(i-1)Δθ, where θ0 is the initial angle of the DOA plane. The steering matrix and reflection vector in (25) are extended to obtain the M×N extended steering matrixΦ and the N×1 extended reflection vector x, which are defined as Φ=[a(θ0)a(θ0+Δθ)⋯a(θ0+(i-1)Δθ)a(θ0+(N-1)Δθ)] and x=s(t)[β(θ0)β(θ0+Δθ)⋯β(θ0+(i-1)Δθ)β(θ0+(N-1)Δθ)]T. Since small number of grids are occupied by the targets, x is a sparse vector with the ith element defined as x(i)=s(t)β(θ0+(i-1)Δθ) if the ith grid is occupied by the target; otherwise, x(i)=0. As a result, the M×1 received complex vector of array observations y could be written as
(26)y=Φx+e,
where e is an M×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 k-term approximation error, and some achievable bounds for the l1 error of the best k-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 is20. 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 (θ0) and angle interval (Δθ) equaling 0° and 3°, respectively. A maximum of L=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 use512 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 {βk=1,k=1,…,4}.The average reconstruction error is adopted to evaluate the reconstruction performance of the methods, which is defined as(27)χaverage=∑i=1NMCχiNMC,
where χi denotes the reconstruction error at the ith Monte Carlo simulation, which is defined as
(28)χi=∥xestimatei-xi∥2∥xi∥2,
where xi and xestimatei represent the true and estimated signal representing the sparse radar scene at the ith 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 Table1 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.
Average reconstruction error
Average RMSE of DOA angles (degree)
SMP
0.04
0.05
OMP
0.2
0.19
ROMP
0.24
0.23
StOMP
0.23
0.25
SP/CoSaMP
0.15
0.13
BAOMP
0.17
0.16
BPDN
0.09
0.07
Capon
0.38
1.1
APES
0.35
0.54
GLRT
0.27
0.17


## 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 ofM antennas, which are used to emit the transmitted signal s(t). The M×1 received complex vector of array observations is defined as f(t)=[f1(t),…,fM(t)]T. 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)f(t)=β(θ)s(t)a(θ)+n(t),
where β(θ) is the reflection coefficient of the hypothetical target, and n(t) is an M×1 complex Gaussian noise vector. a(θ) is the M×1 steering vector, which is defined as
(23)a(θ)=[1ej(2πdsinθ/λ),…,ej(M-1)(2πdsinθ/λ)]T,
where d is the distance between the elements of the arrays and λ denotes wavelength.AssumingD targets are observed with reflection coefficients {βi}i=1D and DOA angles {θi}i=1D, the M×1 received complex vector of array observations can be written as
(24)α(t)=∑i=1Dβ(θi)s(t)a(θi)+γ(t),
where γ(t) is an M×1 complex Gaussian noise vector. Equation (24) could be rewritten as
(25)α(t)=A(θ)S(t,θ)γ(t),
where A(θ)=[a(θ1)a(θ2)⋯a(θD)] is an M×D steering matrix, and S(t,θ)=s(t)[β(θ1)β(θ2)⋯β(θD)]T denotes a D×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 intoN fine grids, each cell generally with the same size Δθ. The ith grid represents the DOA angle of θ0+(i-1)Δθ, where θ0 is the initial angle of the DOA plane. The steering matrix and reflection vector in (25) are extended to obtain the M×N extended steering matrixΦ and the N×1 extended reflection vector x, which are defined as Φ=[a(θ0)a(θ0+Δθ)⋯a(θ0+(i-1)Δθ)a(θ0+(N-1)Δθ)] and x=s(t)[β(θ0)β(θ0+Δθ)⋯β(θ0+(i-1)Δθ)β(θ0+(N-1)Δθ)]T. Since small number of grids are occupied by the targets, x is a sparse vector with the ith element defined as x(i)=s(t)β(θ0+(i-1)Δθ) if the ith grid is occupied by the target; otherwise, x(i)=0. As a result, the M×1 received complex vector of array observations y could be written as
(26)y=Φx+e,
where e is an M×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 k-term approximation error, and some achievable bounds for the l1 error of the best k-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 is20. 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 (θ0) and angle interval (Δθ) equaling 0° and 3°, respectively. A maximum of L=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 use512 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 {βk=1,k=1,…,4}.The average reconstruction error is adopted to evaluate the reconstruction performance of the methods, which is defined as(27)χaverage=∑i=1NMCχiNMC,
where χi denotes the reconstruction error at the ith Monte Carlo simulation, which is defined as
(28)χi=∥xestimatei-xi∥2∥xi∥2,
where xi and xestimatei represent the true and estimated signal representing the sparse radar scene at the ith 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 Table1 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.
Average reconstruction error
Average RMSE of DOA angles (degree)
SMP
0.04
0.05
OMP
0.2
0.19
ROMP
0.24
0.23
StOMP
0.23
0.25
SP/CoSaMP
0.15
0.13
BAOMP
0.17
0.16
BPDN
0.09
0.07
Capon
0.38
1.1
APES
0.35
0.54
GLRT
0.27
0.17


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



---

*Source: 804640-2013-05-12.xml*