Journal of Applied Mathematics

Volume 2013 (2013), Article ID 804640, 8 pages

http://dx.doi.org/10.1155/2013/804640

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

^{1}The School of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China^{2}The MOE Key Lab for Intelligent and Networked Systems, The School of Electronic and Information Engineering, Xian Jiaotong University, Xi’an, Shaanxi 710049, China

Received 14 March 2013; Accepted 16 April 2013

Academic Editor: Xianxia Zhang

Copyright © 2013 Liu Jing et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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 . The resulting adaptation of (2), known as basis pursuit (BP) [1], is formally defined as

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 measurement noise vector. We simply change the constraint on the solution to 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 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 largest components at each iteration. However, in practice some true indices may correspond to a set of small components. The proposed method selects () 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 largest components at each iteration. However, in practice some true indices may correspond to a set of small components. The proposed SMP selects () 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 () 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 .

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

*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 iteration, go through the following steps.(1)*Splitting*. Locate the largest components of the correlation vector at the iteration, and divide them into split sets as
where denotes the split set at the iteration, and denotes the correlation vector at the iteration. Furthermore, denotes the largest magnitude entry to the largest magnitude entry of .(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
where denotes the merged split set at the iteration, and denotes the estimated support set at the 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 () and sets other entries as zero, resulting in a sequence of split estimates as
where denotes the split estimate of the target signal at the 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 , for example,
where denotes the pruned split set at the iteration. (5)*Split Sets Merging*. The pruned split sets are merged to form a final merged set, , as
(6)*Estimation*. An interim estimate of the original signal is calculated based on the final merged set using least-square algorithm as
where denotes the interim estimate of the original signal at the 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
where denotes the estimated support set at the iteration. (8)*Estimation*. Obtain the final estimate at each iteration, based on using least-square algorithm as
where denotes the final estimate at the 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:
where denotes the residual at the iteration. (10) If , perform the correlation calculation , and then go to step of the iteration; otherwise, set 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
**then the signal approximation is -sparse and**In particular,
**
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
**
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 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 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 guarantees if .*

*Proof. *Considering the monotonicity of : for any two integers , according to Lemma 1 in [7]. So we have provided that .

Proposition 5. *Step 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 , resulting in . Since is set large enough to expand the searching area, is set as in the proposed method. We then have according to . Propositions 3~5 are based on Theorem 2, which has a rigid setting for . We can relax the range of to . As a result, we need not know the exact value of and can select an integer randomly from 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 , which could be approximated as when is chosen as a small value (e.g. 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
**
where .*

Proposition 7. *The number of iterations of the SMP algorithm is upper bounded by
*

*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 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 given the length of the original signal () and the length of the measurement vector (). 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 , , , , as the input parameters.

The signal sparsity level is varied from to . 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 , that is, . The frequency of exact reconstruction () is used to evaluate the reconstruction performance of the different methods, which is defined as where denotes the number of exact reconstructions for each algorithm given a fixed , and 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 , 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.

##### 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 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 where is the reflection coefficient of the hypothetical target, and is an complex Gaussian noise vector. is the steering vector, which is defined as where is the distance between the elements of the arrays and denotes wavelength.

Assuming targets are observed with reflection coefficients and DOA angles , the received complex vector of array observations can be written as where is an complex Gaussian noise vector. Equation (24) could be rewritten as where is an steering matrix, and denotes a 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 grid represents the DOA angle of , 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 extended reflection vector , which are defined as and . Since small number of grids are occupied by the targets, is a sparse vector with the element defined as if the grid is occupied by the target; otherwise, . As a result, the received complex vector of array observations could be written as where is an 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 . The antennas transmit independent orthogonal quadrature phase shift keyed (QPSK) waveforms and the carrier frequency is GHz. The SNR of the measurement noise is set to a fixed value (20 dB). The range of the DOA plane is , which is divided into cells with the initial angle () and angle interval () equaling and , respectively. A maximum of 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 snapshots each. Five hundred Monte Carlo simulations are carried out, and in each trial four targets locate randomly within the DOA range of , and the corresponding reflection coefficients are set as .

The average reconstruction error is adopted to evaluate the reconstruction performance of the methods, which is defined as where denotes the reconstruction error at the th Monte Carlo simulation, which is defined as where and represent the true and estimated signal representing the sparse radar scene at the 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 snapshots.

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

#### References

- S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,”
*SIAM Journal on Scientific Computing*, vol. 20, no. 1, pp. 33–61, 1998. View at Publisher · View at Google Scholar · View at MathSciNet - S. Boyd and L. Vandenberghe,
*Convex Optimization*, Cambridge University Press, Cambridge, UK, 2004. View at MathSciNet - J. A. Tropp and S. J. Wright, “Computational methods for sparse solution of linear inverse problems,”
*Proceedings of the IEEE*, vol. 98, no. 6, pp. 948–958, 2010. View at Publisher · View at Google Scholar · View at Scopus - S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,”
*IEEE Transactions on Signal Processing*, vol. 41, no. 12, pp. 3397–3415, 1993. View at Publisher · View at Google Scholar · View at Scopus - J. A. Tropp, “Greed is good: algorithmic results for sparse approximation,”
*IEEE Transactions on Information Theory*, vol. 50, no. 10, pp. 2231–2242, 2004. View at Publisher · View at Google Scholar · View at MathSciNet - D. Needell and R. Vershynin, “Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit,”
*Foundations of Computational Mathematics*, vol. 9, no. 3, pp. 317–334, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - W. Dai and O. Milenkovic, “Subspace pursuit for compressive sensing signal reconstruction,”
*IEEE Transactions on Information Theory*, vol. 55, no. 5, pp. 2230–2249, 2009. View at Publisher · View at Google Scholar · View at MathSciNet - D. Needell and J. A. Tropp, “CoSaMP: iterative signal recovery from incomplete and inaccurate samples,”
*Applied and Computational Harmonic Analysis*, vol. 26, no. 3, pp. 301–321, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - H. Huang and A. Makur, “Backtracking-based matching pursuit method for sparse signal reconstruction,”
*IEEE Signal Processing Letters*, vol. 18, no. 7, pp. 391–394, 2011. View at Publisher · View at Google Scholar · View at Scopus - J. Capon, “High resolution frequency-wave number spectrum analysis,”
*Proceedings of the IEEE*, vol. 57, no. 8, pp. 1408–1418, 1969. View at Google Scholar · View at Scopus - J. Li and P. Stoica, “An adaptive filtering approach to spectral estimation and SAR imaging,”
*IEEE Transactions on Signal Processing*, vol. 44, no. 6, pp. 1469–1484, 1996. View at Google Scholar · View at Scopus - E. J. Kelly, “An adaptive detection algorithm,”
*IEEE Transactions on Aerospace and Electronic Systems*, vol. 22, no. 2, pp. 115–127, 1986. View at Google Scholar · View at Scopus - SparseLab, http://dsp.rice.edu/cs.
- A. C. Gürbüz, J. H. McClellan, and V. Cevher, “A compressive beamforming method,” in
*Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP '08)*, pp. 2617–2620, Las Vegas, Nev, USA, April 2008. View at Publisher · View at Google Scholar · View at Scopus - I. Stojanovic, W. C. Karl, and M. Cetin, “Compressed sensing of mono-static and multi-static SAR,” in
*Algorithms for Synthetic Aperture Radar Imagery XVI*, Proceedings of SPIE, Orlando, Fla, USA, April 2009. View at Publisher · View at Google Scholar · View at Scopus - J. H. G. Ender, “On compressive sensing applied to radar,”
*Signal Processing*, vol. 90, no. 5, pp. 1402–1414, 2010. View at Publisher · View at Google Scholar · View at Scopus - Y. Yu, A. P. Petropulu, and H. V. Poor, “Measurement matrix design for compressive sensing-based MIMO radar,”
*IEEE Transactions on Signal Processing*, vol. 59, no. 11, pp. 5338–5352, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - J. Zhang, D. Zhu, and G. Zhang, “Adaptive compressed sensing radar oriented toward cognitive detection in dynamic sparse target scene,”
*IEEE Transactions on Signal Processing*, vol. 60, no. 4, pp. 1718–1729, 2012. View at Publisher · View at Google Scholar · View at MathSciNet - Y. Chi, L. L. Scharf, A. Pezeshki, and A. R. Calderbank, “Sensitivity to basis mismatch in compressed sensing,”
*IEEE Transactions on Signal Processing*, vol. 59, no. 5, pp. 2182–2195, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - Y. Bar-Shalom and T. E. Fortmann,
*Tracking and Data Association*, Academic Press, Orlando, Fla, USA, 1988. View at MathSciNet