- Email: [email protected]
A spatial-domain multiresolutional particle filter with thresholded wavelets
ARTICLE IN PRESS
Signal Processing 87 (2007) 1384–1401 www.elsevier.com/locate/sigpro
A spatial-domain multiresolutional particle filter with thresholded wavelets Lang Honga,, Devert Wickerb a
Department of Electrical Engineering, Wright State University, Dayton, OH 45435, USA Target Recognition Branch, Air Force Research Laboratory, SNAT, WPAFB, OH 45433, USA
b
Received 10 July 2006; received in revised form 19 July 2006; accepted 5 December 2006 Available online 10 January 2007
Abstract Particle filters have been proven to be very effective for nonlinear/non-Gaussian filtering. However, the most notorious disadvantage of a particle filter is its formidable computational complexity, since hundreds (even thousands) of particles are usually needed to achieve a required approximation accuracy. It has also been proven that one of the techniques of truly solving a computational problem is multiresolutional processing, both in temporal and spatial domains. Therefore, in this paper we propose a multiresolutional particle filter in the spatial domain using thresholded wavelets to reduce significantly the number of particles, meanwhile maintaining the full strength of a particle filter. r 2007 Elsevier B.V. All rights reserved. Keywords: Particle filters; Nonlinear/non-Gaussian estimation; Computational efficiency; Multiresolutional processing techniques; Wavelets
1. Introduction Particle filters, falling into the category of sequential Monte Carlo methods, are effective in nonlinear and non-Gaussian filtering. A particle filter estimates a posterior density function pðxk jZ k Þ by an empirical approach ^ k jZ k Þ ¼ pðx
N 1X dx ðxkðiÞ Þ, N i¼1
(1)
where xkðiÞ are iid (independently and identically distributed) random samples of pðxk jZ k Þ and Zk represents the accumulative measurements. The Corresponding author.
E-mail addresses: [email protected] (L. Hong), [email protected] (D. Wicker).
associated integration can also be achieved empirically as Z N 1X f ðxðiÞ Þ. (2) Iðf t Þ ¼ f t ðxk Þpðxk jZ k Þ dxk N i¼1 t k Although operations in Eqs. (1) and (2) are straightforward and can be applied to any nonlinear and non-Gaussian density, there is a problem: sampling directly from a posterior function is not always possible. Researchers have found a solution: instead of sampling the posterior function, samples can be drawn from a proposal function whose support includes that of the posterior function. This method is called importance sampling (IS) or sequential importance sampling (SIS) when the weight updates are performed in a recursive manner. However, IS/SIS is not problem free. As time elapses, the
0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2006.12.005
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
distribution of importance weights becomes more and more skewed. After a few time steps, only a few particles have non-zero weights. Another fix has also been proposed: a selection step called importance resampling is introduced to eliminate the particles with low weights and multiply the particles having high weights. More recently, many improvements have been made and the technology of particle filters is becoming mature, which includes bootstrap particle filters, extended-Kalman particle filters, unscented particle filters, etc. [1,2]. Ref. [3] presents an excellent collection of articles on particle filters. When applied to target tracking, both references [4,5] give a good comparison of different nonlinear filtering approaches, including particle filters. Although particle filters are effective for nonlinear/non-Gaussian filtering, the most notorious disadvantage of a particle filter is its formidable computational complexity, since hundreds (even thousands) of particles are usually needed to achieve a decent density approximation accuracy. Little effort has been attempted with limited success on improving the computational efficiency of a particle filter [4,6,7]. It is well known that there are two fundamental approaches in achieving computational efficiency in engineering applications: parallel processing and multiresolutional processing. Parallel processing divides a computational task into many subtasks which can then be executed simultaneously with multiple processors, while multiresolutional processing prioritizes the computational resources according to the need to maintain the ‘‘details’’ of signals. Many engineering applications cannot be implemented in parallel due to the nature of a single computational thread. On the other hand, multiresolutional processing has a broader application range. This paper presents a computationally efficient particle filter using a wavelet-based multiresolutional estimation technique developed over the last decade by the first author [8–10]. In general, multiresolutional processing can be applied to a particle filter in either temporal or spatial domains. A multiresolutional particle filter is called a multirate particle filter if it is derived by applying multiresolutional processing in a temporal domain. On the other hand, applying multiresolutional processing in a spatial domain in designing a particle filter results in a multires particle filter. In our previous work, we have demonstrated significant computational savings by using a multirate multiple model particle filter [11]. This paper is focused on computational savings with a multires
1385
particle filter. A particle filter approximates a density function by a set of particles with weights and a multires particle filter further decomposes the weights into lowpassed and highpassed components. By thresholding the highpassed particle weights, a much smaller number of particles can represent the density function with a comparable approximation accuracy. A density can be approximated by a set of densities with different resolutions (smoothness), where a smoother density requires a smaller number of particles. As a result, the total number of particles could be much smaller than a regular particle filter (we call this a unires particle filter), while maintaining a comparable approximation accuracy. The rest of the paper is organized as follows. Section 2 formulates a multires particle filter and Section 3 presents detailed implementations of a multires particle filter with simulation examples. Section 4 concludes the paper. 2. Formulation of a multires particle filter For a nonlinear filtering problem, the system and measurement models are given by xk ¼ f ðxk1 Þ þ vk1 ,
ð3Þ
zk ¼ hðxk Þ þ wk ,
ð4Þ
where f ðÞ and hðÞ are system dynamic and measurement functions. Variables vk1 and wk are governed by the following known densities, pv and pw : pðxk jxk1 Þ ¼ pv ðxk f ðxk1 ÞÞ, pðzk jxk Þ ¼ pw ðzk hðxk ÞÞ.
ð5Þ ð6Þ
The initial state density, px ðx0 Þ, is also assumed known which can be decomposed into a multires structure with j 1 levels X XX px ðx0 Þ ¼ aj1 ;k fj1 ;k ðx0 Þ þ bj;k cj;k ðx0 Þ (7) k
jXj 1
k
with scaling and wavelet coefficients given by [12] aj1 ;k ¼ hpx ðx0 Þ; fj 1 ;k ðx0 Þi,
ð8Þ
bj;k ¼ hpx ðx0 Þ; cj;k ðx0 Þi,
ð9Þ
where ffj1 ;k ¼ 2j 1 =2 fð2j1 x kÞ; k 2 Zg, fcj;k ¼ 2j=2 cð2j x kÞ; k 2 Z; jXj 1 g, and fðxÞ and cðxÞ are the scaling and mother wavelet functions, respectively. In the form of
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1386
particle representation, we have px ðx0 Þ ¼
N0 X
oix0 dðx0 xi0 Þ,
(10)
i¼1
where particle xi0 is sampled from px ðx0 Þ, i.e., xi0 px ðx0 Þ and oix0 is the corresponding weight. The inner products defined by Eqs. (8) and (9) that require nonlinear integration become a^ j1 ;k ¼
N0 X
ol 1 ¼
N 0 =2 X
i
(18)
i
(19)
oil 1 dðx0 x01 Þ
i¼1
and oix0 fj1 ;k ðxi0 Þdðx0
xi0 Þ
(11)
i¼1
oh1 ¼
and b^ j;k ¼
where only level 1 decomposition is presented. Not considering the edge effect, the process of convolution and downsampling results in N20 non-zero coefficients which can be described in a summation form
N 0 =2 X
oil 1 dðx0 x01 Þ,
i¼1 N0 X
i
oix0 cj;k ðxi0 Þdðx0
xi0 Þ.
(12)
i¼1
In this paper, we will use Harr wavelets as a vehicle for demonstration.1
where x01 is the result of lowpassed and downsampled xi0 . Level 2 decomposition can be derived from ol 1 by ol 2 ½i ¼
þ1 X
ol 1 ½i h½2i k
(20)
ol 1 ½i g½2i k,
(21)
k¼1
2.1. Harr wavelets
and
The scaling function is a simple rectangular function 1 if 0ptp1; fðtÞ ¼ (13) 0 otherwise with only two non-zero coefficients hð0Þ ¼ hð1Þ ¼ p1ffiffi2 and the wavelet function is 8 if 0pto0:5; > <1 (14) cðtÞ ¼ 1 if 0:5ptp1; > : 0 otherwise with two non-zero coefficients gð0Þ ¼ p1ffiffi2 and gð1Þ ¼ p1ffiffi2. A filterbank is used for the discrete wavelet transform (DWT). Rewrite Eq. (10) as a signal sequence px ðx0 Þ ¼ fox0 ½ig;
i ¼ 1; . . . ; N 0 .
(15)
Eqs. (11) and (12) can be expressed as a discrete convolution with a downsampling process ol 1 ½i ¼ a^ 1 ½i ¼
þ1 X
ox0 ½i h½2i k
(16)
oh2 ½i ¼
þ1 X k¼1
where only N40 non-zero coefficients are obtained. Other levels of decomposition can be obtained similarly. The major difference from the traditional wavelet transform in Eqs. (11) and (12) is that xi0 ; i ¼ 1; . . . ; N are random samples, not a predefined and evenly spaced sequence. To overcome the randomness, we sample all N 0 particles first; then filtering is applied, but the samples are still unevenly spaced. Because of the downsampling operation, the numbers of transformed particles ol j1 and ohj are drastically reduced with the number of levels. With the scaling and wavelet coefficients estimated by Eqs. (16) and (17), the density (weights) reconstruction can be performed by ox0 ½i ¼
þ1 X
ol 1 ½kh0 ½2k i þ oh1 ½kg0 ½2k i.
k¼1
k¼1
(22)
and oh1 ½i ¼ b^ 1 ½i ¼
þ1 X
ox0 ½i g½2i k,
(17)
k¼1 1
The multires technology developed here can be applied to any wavelet.
A simple thresholding can be carried out on the wavelet coefficients to reduce the number of coefficients. Let ~ h1 ½i ¼ oh1 ½i if oh1 ½iXT s ; o
~ h1 ½i ¼ 0. otherwise o
(23)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
A thresholded density estimation can then be achieved as þ1 X
ox0 ½i ¼
~ h1 ½kg0 ½2k i. ol 1 ½kh0 ½2k i þ o
k¼1
(24) T s is an application dependent parameter which needs to be adjusted for each application depending on available computational capability. The larger T s is, the fewer number of surviving coefficients but larger approximation errors. Example 1. Assume density px ðx0 Þ is given by the following bi-modal function: px ðx0 Þ ¼ 0:6Nð10; 22 Þ þ 0:4Nð20; 42 Þ,
(25)
where NðÞ denotes a normal distribution function. px ðx0 Þ is approximated by Eq. (10) whose particle weights are plotted in Fig. 1 with the number of particles being 1000. A three-level multires decomposition for the particle weights is performed and the results are shown in Fig. 2. The original particle weights are now represented by multiresolutional weights in Fig. 2, where (a) and (b) are level-1 lowpass-filtered and highpass-filtered weights, (c) and (d) are level-2 lowpass-filtered and highpass-filtered weights, and (e) and (f) are level-3 lowpass-filtered and highpass-filtered weights. In the reconstruction process, weights from (b), (d)–(f) completely reconstruct the weights in Fig. 1. When a threshold is applied to the highpassed particle weights, the numbers of particles can be drastically reduced at each level without significantly degrading reconstruction accuracy. Figs. 3 and 4 show the thresholding results, where Fig. 3 shows that only a small number of highpass-filtered particle weights
are used in the reconstruction with the threshold being 0.001. The results for threshold 0.002 are given in Figs. 5 and 6. Thus far, the decomposition is performed for the sequence where the independent variable is the sequence index i. If the independent variable is a continuous variable, such as x, the decomposition is handled differently. Assume density px can be approximated by px ¼
N X
oi dðx xi Þ,
(26)
i¼1
where xi is sampled from px , i.e., xi px , and P N i¼1 oi ¼ 1. Density px is a continuous function of x, not i. To better handle Eq. (26) in multiresolution decomposition, we can rewrite it in a vector form px ¼ dT o ,
(27)
where d ¼ ½dðx x1 Þ; . . . ; dðx xN ÞT , o ¼ ½o1 ; . . . ; oN T .
ð28Þ
Let o ¼ TT ½ol j ; ohj ; . . . ; ohj 1
1
1 þL
px ¼ dT TT olh ¼ ðT d ÞT olh ¼ dTlh olh ,
0.008 0.006 0.004 0.002 0 0
5
10
15
20
25
30
35
x Fig. 1. The weights of a particle approximation of the density defined in Eq. (25).
(29)
(30)
where 1
0.01
T ¼ TT olh
be the multiresolution decomposition of sample weights, where T is an orthogonal wavelet decomposition transform matrix, such that T1 ¼ TT . There are L þ 1 levels in the decomposition. Substituting Eq. (29) into Eq. (27) results in
dlh ¼ T d ¼ ½dl j ; dhj ; . . . ; dhj
0.012
particle weights
1387
1
1 þL
T .
(31)
The multiresolution decomposition in Eq. (30) is applied to both weights and the independent variable.2 There are two methods of implementing Eq. (31): implicit and explicit. The implicit method embeds the wavelet transformation into a complicated variable structure and no inverse transform is needed when performing density reconstruction. On the other hand, the explicit method uses a simple variable structure, but an inverse transform is needed in density reconstruction. 2
To be precise, multiresolution decomposition is applied to the d function of the independent variable.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1388
b
a
x 10-3 4
0.008
level 1 highassed
level 1 lowpassed
0.01
0.006 0.004 0.002 0
2 0 -2 -4
0
5
10
15
20
25
30
35
0
5
10
15
x
c
d
x 10-3
level 2 highpassed
level 2 lowpassed
25
30
35
20
25
30
35
20
25
30
35
x 10-3 4
8 6 4 2 0
2 0 -2 -4
0
5
10
15
20
25
30
35
0
5
10
15
x
e
x
f
x 10-3
1 level 3 highpassed
6 level 3 lowpassed
20 x
4
2
0
x 10-3
0.5 0 -0.5 -1 -1.5
0
5
10
15
20
25
30
35
0
5
10
x
15 x
Fig. 2. A three-level multires decomposition of the particle weights.
Implicit method: To make an equivalent variable for each dl and dh in Eq. (31), we introduce the following variable structure:
and sgn is a sign function ( 1; x 2 first half of fxhj g; sgnðfxhj gÞ ¼ 1; x 2 second half of fxhj g:
dl j ¼ hj 1 dðx fxl gÞ,
(34) With this implicit method, the density reconstruction can be performed as follows: Nl X px ¼ oil dðx fxl gi Þ
1
pffiffiffi 2 j , ð32Þ dhj ¼ h dðx fxhj gÞsgnðfxhj gÞ; jXj 1 ; h ¼ 2 where fxl g and fxhj g are collections of sample points
i N
fxl g ¼ fsamples associated with ol g fxhj g ¼ fsamples associated with ohj g
and
þ ð33Þ
hj XX
jXj 1
i
oihj dðx fxhj gi Þsgnðfxhj gi Þ.
ð35Þ
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
a
b
x 10-3 4
level 1 highassed
level 1 lowpassed
0.01 0.008 0.006 0.004 0.002 0
2 0 -2 -4 -6
0
5
10
15
20
25
30
35
0
5
10
15
x
c
d level 2 highpassed
level 2 lowpassed
6 4 2 0
30
35
20
25
30
35
20
25
30
35
x 10-3 2 1 0 -1 -2
0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
x 10-3 6
x 10-3 2
level 3 highpassed
level 3 lowpassed
25
3
8
4
2
0
20 x
x 10-3
e
1389
0
5
10
15
20
25
30
35
1 0 -1 -2
0
5
10
15
x
x
Fig. 3. A three-level multires decomposition of the particle weights with thresholding where the threshold is 0.001.
Explicit method: It can be seen in Eq. (35) that no explicit inverse transform is needed in density reconstruction, but the variable structure in Eq. (32) may pose a difficulty in implementation. Other than a complicated variable structure, we can introduce transformed variables as 3 2 2 13 xl x 6 x 7 6 x2 7 6 h1 7 6 7 7 6 7 (36) xslh ¼ 6 .. 7 ¼ T6 6 .. 7. 6 . 7 4 5 . 5 4 xhL2L xN
Then by using dlh ¼ T d ¼ Tdðx TT xslh Þ,
(37)
we have px ¼ ðTdðx TT xslh ÞÞT olh T
¼ ðTdðT ðxlh
xslh ÞÞÞT olh
ð38Þ ð39Þ
which contains explicit transformation T. Eq. (38) can be rewritten as px ¼ pT ðdlh ; olh Þ, where pT is a transformation mapping function.
(40)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1390
b
0.01 0.008
original level 2 lowpassed
original level 2 lowpassed
a
0.006 0.004 0.002
0.01 0.008 0.006 0.004 0.002
0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
d
0.01 0.008
original level 1 lowpassed
original level 1 lowpassed
c
0.006 0.004 0.002
30
35
20
25
30
35
20
25
30
35
0.008 0.006 0.004 0.002 0
0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
0.01 0.008
0.01
original level 0 lowpassed
original level 0 lowpassed
25
0.01
0
e
20 x
0.006 0.004 0.002
0.008 0.006 0.004 0.002
0
0 0
5
10
15
20
25
30
35
0
5
x
10
15 x
Fig. 4. A density reconstruction with thresholding where the threshold is 0.001.
If we name each level’s density component3 in Eq. (30) by pjx1l ¼ dTl ol
(41)
for the lowpassed density component at level j 1 for the implicit method or pjx1l ¼ pT ðdl ; ol Þ
(42)
for the explicit method, and pjxh ¼
X
dThjk ohjk ;
j ¼ j1; . . . ; L
(43)
k 3 We call them density components, since they are not regular densities anymore, due to the fact that the integration of each component is not one.
for the highpassed density component at level j using the implicit method, or X pT ðdhjk ; ohjk Þ; j ¼ j 1 ; . . . ; L (44) pjxh ¼ k
using the explicit method, we could rewrite Eq. (30) as X px ¼ pjx1l þ pjxh . (45) jXj 1
With simple thresholding, we can significantly reduce the number of particles without sacrificing the performance. Thresholding can be applied to both the implicit and explicit methods. For the implicit method, thresholding is applied to the density components 8 T < fdhjk ohjk gi ; jfdThjk ohjk gi j4T s ; T i ~ hjk g ¼ fdhjk o (46) : 0; jfdThjk ohjk gi jpT s
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
a
b
x 10-3 4
level 1 highpassed
level 1 lowpassed
0.015 0.01 0.005 0
2 0 -2 -4
0
5
10
15
20
25
30
35
0
5
10
15
x
c
d
30
35
20
25
30
35
20
25
30
35
5 level 2 highpassed
level 2 lowpassed
25
x 10-3
8 6 4 2
4 3 2 1
0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
x 10-3
1 level 3 highpassed
8 level 3 lowpassed
20 x
x 10-3
e
1391
6 4 2
0.5 0
-0.5
0
-1 0
5
10
15
20
25
30
35
0
5
10
15
x
x
Fig. 5. A three-level multires decomposition of the particle weights with thresholding where the threshold is 0.002.
and for the explicit method, thresholding is directly applied to the wavelet coefficients ( ohjk ; johjk j4T s ; ~ hjk ¼ o (47) 0; johjk jpT s : Thus, we have X ~ hjk ; dThjk o pjxh ¼
Example 2. Let the wavelet transform be 3 2 1=2 1=2 1=2 1=2 7 6 1=2 1=2 7 6 1=2 1=2 pffiffiffi pffiffiffi 7, T¼6 6 2=2 2=2 0 0 7 4 pffiffiffi pffiffiffi 5 0 0 2=2 2=2
(50)
and a 4-point density approximation be j ¼ j1 ; . . . ; L
(48)
k
for the implicit method, and X ~ hjk Þ; j ¼ j 1 ; . . . ; L pjxh ¼ pT ðdhjk ; o k
for the explicit method.
½x1 x2 x2 x4 ¼ ½2 10 4 6T , ½o1 o2 o3 o4 ¼ ½16
(49)
1 5 1T 4 12 6
ð51Þ
and 5 px ¼ 16 dðx 2Þ þ 12 dðx 4Þ þ 16 dðx 6Þ þ 14 dðx 10Þ.
(52)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1392
b
0.01
0.01
original level 2 lowpassed
original level 2 lowpassed
a
0.008 0.006 0.004
0.008 0.006 0.004 0.002
0.002 0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
d
0.01
original level 1 lowpassed
original level 1 lowpassed
c
0.008 0.006 0.004
30
35
20
25
30
35
20
25
30
35
0.008 0.006 0.004 0.002
0.002
0 0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
0.01
0.01
original level 0 lowpassed
original level 0 lowpassed
25
0.01
0
e
20 x
0.008 0.006 0.004
0.008 0.006 0.004 0.002
0.002
0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
x
Fig. 6. A density reconstruction with thresholding where the threshold is 0.002.
Implicit method: From Eqs. (29) and (31), we have 2
dl
3
7 6 6 dh 1 7 7 6 7 ¼ Td 6 6 dh21 7 5 4 dh22 2 3 1=2½dðx 2Þ þ dðx 4Þ þ dðx 6Þ 6 7 6 7 þdðx 10Þ 6 7 6 7 6 1=2½dðx 2Þ dðx 4Þ þ dðx 6Þ 7 7 6 ¼6 7 7 6 þdðx 10Þ 7 6 pffiffiffi 7 6 7 6 2 =2½dðx 2Þ þ dðx 4Þ 5 4 pffiffiffi 2=2½dðx 6Þ þ dðx 10Þ ð53Þ
and 2 1 3 3 2 ol 2 6 17 6o 7 6 12 7 6 h1 7 6 pffiffi 7 7 6 7. 6 oh 7 ¼ T o ¼ 6 6 82 7 4 21 5 4 pffiffi 5 2 oh22
(54)
24
Then px ¼ dl ol þ dh1 oh1 þ dh21 oh21 þ dh22 oh22 ð55Þ 1 ¼ 4 ½dðx 2Þ þ dðx 4Þ þ dðx 6Þ þ dðx 10Þ 1 24 ½dðx 2Þ dðx 4Þ þ dðx 6Þ þ dðx 10Þ þ 18 ½dðx 2Þ 1 þ dðx 4Þ þ 24 ½dðx 6Þ þ dðx 10Þ. ð56Þ
For this example, fxl g ¼ f2; 4; 6; 10g, fxh1 g ¼ f2; 4; 6; 10g, fxh21 g ¼ f2; 4g and fxh22 g ¼ f6; 10g. If we set
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1393 j
1 the threshold as 24 , applying thresholding to Eq. (56) results in
Figs. 7 and 8 show the reconstruction of px from px1l and pjxh with and without thresholding using the implicit method.
px dl ol þ dh1 oh1 þ dh21 oh21 þ dh22 oh22
Explicit method: Using Eq. (36), we can find the following transformed variables: pffiffiffi pffiffiffi xlh ¼ ½11 5 2 2 2T . (59)
¼
1 4 ½dðx 2Þ þ dðx 4Þ þ dðx þ 18 ½dðx 2Þ þ dðx 4Þ.
ð57Þ
6Þ þ dðx 10Þ ð58Þ
5 5 12 12
px
Threshold = 0 4 multires particles reconstruct original 4 particles
1 1 6 6
1 1 4 4
1 1 6 6
x
px
h2
2 1 24
6 4
- 1 8
6
8
10
4 particles
1 8 2
4
8 - 1 24 2 particles
Direct sum; no inverse transformation is needed; but complicated variable structure is needed
+
x
10
px
px
l
h1
1 4 1 24 2
6 4
8
10
x
x
- 1 24
4
2
6
1 particle
8
10
1 particle j
Fig. 7. Multires density reconstruction from px1l and pjxh using the implicit method without thresholding.
px Threshold = 1 24 2 multires particles approximate original 4 particles
1 6 1 8
5 12 3 8
x
px
2
h2
2
- 1 8
8
10
+ Direct sum; no inverse transformation is needed; but complicated variable structure is needed
x 6 1 particle
8 6 4 particles
4
1 8
4
1 1 4 4
1 1 4 6
10
px
px
l
h1
1 4 x 2
4
6 0 particle
8
x
10
2 j
4
6 1 particle
8
10
Fig. 8. Multires density reconstruction from px1l and pjxh using the implicit method with thresholding.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1394
The reconstruction of the density can be performed separately by variable and weight reconstructions, which are
a density in a multires decomposition form X pjxh , pðxk1 jz1:k1 Þ ¼ pjx1l þ
x ¼ ½2 4 6 10T ¼ TT xlh
where
(60)
and o¼
1
1 T
5 1 6 12 6 4
pjx1l
T
¼ T olh .
(61)
The density reconstruction is then px ¼
1 6 dðx
2Þ þ
k1
5 12 dðx
4Þ þ
1 6 dðx
6Þ þ
1 4 dðx
10Þ
whichpisffiffi illustrated in Fig. 9. If the threshold is set as T s ¼ 242, then we have pffiffi T 2 1 olh ¼ ½12 12 (63) 8 0 ,
¼ dTl k1 ol k1
and pjxh
(62)
¼ k1
X m
dThjm
(66)
k1
~ hjm o
k1
;
j ¼ j 1 ; . . . ; L.
(67)
In the following, we will derive a complete update cycle for multires densities. 3.1. The implicit method Multires propagation:
which results in o ¼ ½16
k1
(65)
k1
jXj 1
5 5 5 T 12 24 24 .
(64)
Figs. 9 and 10 demonstrate the results of density reconstruction using the explicit method, with and without thresholding. 3. Multires bootstrap particle filter algorithm
pðxk jz1:k1 Þ Z ¼ pðxk jxk1 Þpðxk1 jz1:k1 Þ dxk1 Z ¼ pðxk jxk1 Þpjx1l dxk1 Z X þ pðxk jxk1 Þ pjxh dxk1 ¼
Nl Z X
pðxk jfxl gik1 Þoil k1 dðxk1 fxl gik1 Þ dxk1
i¼1
5 5 12 12
px
1 1 4 4
1 1 6 6
1 1 6 6
Threshold = 0 4 multires particles reconstruct original 4 particles
x
ωh2
24
8
6
10
4 particles 2 24
TT
xh2
2 22 2 4
6
8
Need an inverse transform
10
2 particles 1 2
ωl
ωh1
xh1
5 2
4 6 8 1 1 particle 12
10
ð69Þ
jXj 1
This section presents a major development of a spatial domain multiresolution (multires) particle filter for computationally efficient processing. Given
2 8
ð68Þ
xl 2
4
6
8
10 11
1 particle
Fig. 9. Multires density reconstruction from olh and xlh using the explicit method without thresholding.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
5 5 12 12
px 2 Threshold = 24 3 multires particles reconstruct original 4 particles
1395
1 6
1 4
1 5 6 24
1 6
5 24 x
ωh2 2
4
8
6
10
4 particles
2 8
TT
xh2 4
2 2
6 8 1 particle
Need an inverse transform
10
1 2
ωl
ωh1
xh1
5 2
4 6 8 - 1 12 1 particle
10
xl 2
4
6 1 particle
8
10 11
Fig. 10. Multires density reconstruction from olh and xlh using the explicit method with thresholding.
hj XZ XX
N
þ
jXj 1 i¼1
m
The propagated state estimate can be derived as pðxk jfxhjm gik1 Þoihk1
dðxk1 fxhjm gik1 Þsgnðfxhjm gik1 Þ dxk1 ¼
Nl X
x^ kjk1 ¼ ð70Þ
Nl X
pðxk jfxl gik1 Þoil k1
þ
jXj 1 i¼1
m
~ ihk1 pðxk jfxhjm gik1 Þo
sgnðfxhjm gik1 Þ ¼
Nl X
oil k1 dðxk
ð71Þ
2=3
jXj 1 i¼1
m
vk Nð0; s2v Þ, 2=3
~ ihk1 o
sgnðfxhjm gikjk1 Þ,
dðxk
ð73Þ
(74)
where sv ¼ 1. The propagation density is pðxk jxk1 Þ ¼ Nðxk xk1 ; s2v Þ.
N
þ
~ ihk1 fxhjm gikjk1 o
Example 3. Continuing with Example 2. Assume the system model is xkþ1 ¼ xk þ vk ;
fxl gikjk1 Þ
i¼1 hj X XX
m
sgnðfxhjm gikjk1 Þ.
N
þ
N hj jX 1 þL X X j¼j 1 i¼1
i¼1 hj X XX
oil k1 fxl gikjk1
i¼1
fxhjm gikjk1 Þ ð72Þ
where fxl gikjk1 pðxk jfxl gik1 Þ and fxhjm gikjk1 pðxk j fxhjm gik1 Þ. Since both fxl gikjk1 and fxhjm gikjk1 are associated with a same set of particles at the original resolution level, only on one propagated set (either fxl gikjk1 or fxhjm gikjk1 ) is obtained from sampling and others can be derived from the sample set.
(75)
Take fxl gk1 ¼ f2; 4; 6; 10g for propagation. By sampling pðxk jfxl gik1 Þ, we have fxl gkjk1 ¼ f1:5; 2:8; 3:5; 5:4g. Due to the relationship defined by T in Eq. (50), we can derive other propagated sets fxh1 gkjk1 ¼ f1:5; 2:8; 3:5; 5:4g;
and
fxh22 gkjk1 ¼ f3:5; 5:4g.
fxh21 gkjk1 ¼ f1:5; 2:8g
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1396
Multires update: The update of propagated multires density is given by
There will be a selection process at each level if necessary.
pðxk jz1:k Þ Z 1 ¼ pðzk jxk Þ pðxk jxk1 Þpðxk1 jz1:k1 Þ dxk1 C
3.2. The explicit method
1 ¼ pðzk jxk Þpðxk jz1:k1 Þ C " N l 1 X pðzk jxk Þoil k1 dðxk fxl gikjk1 Þ ¼ C i¼1
ð76Þ
Explicit multires propagation: The unires particle propagation is given by
ð77Þ
xikþ1jk ¼ f ðxik Þ þ vik ;
XXX jXj 1 i¼1
m
pðxk jz1:k1 Þ ¼
~ ihk1 dðxk fxhjm gikjk1 Þ pðzk jxk Þo
¼
N X
oik1 dðxk xikjk1 Þ
i¼1 dTkjk1 ok1
#
¼
ð78Þ
ð86Þ
X
ðdjkjk1 ÞT ojk1
ð87Þ
j¼1
" N l 1 X pðzk jfxl gikjk1 Þoil k1 ¼ C i¼1
¼
L N=2 X
j¼1
ðTdðTT ðxjlhk xjlhkjk1 ÞÞÞT ojlhk1 , ð88Þ
N
þ
jXj 1 i¼1
m
~ ihk1 pðzk jfxhjm gikjk1 Þo
¼
where 2
#
sgnðfxhjm gikjk1 Þ Nl X
ð79Þ
oil k dðxk fxl gikjk1 Þ
i¼1 N hj
þ
XXX jXj 1 i¼1
m
~ ihk dðxk fxhjm gikjk1 Þ o
sgnðfxhjm gikjk1 Þ,
ð80Þ
where the updated weights are oil k ¼
1 pðzk jfxl gikjk1 Þoil k1 C
(81)
and ~ ihk o
(82)
The updated estimate vector can be derived from
x^ kjk ¼
oil k fxl gikjk1
þ
i¼1
sgnðfxhjm gikjk1 Þ.
N hj jX 1 þL X X j¼j 1 i¼1
6 6 xjlhkjk1 ¼ T6 6 4
f ðx1k1 Þ
3j
2
v1k1
3j
7 6 7 7 6 7 7 þ T6 ... 7 7 6 7 5 4 L 5 2L 2 f ðxk1 Þ vk1 .. .
ð89Þ
¼ T f ðTT xjlhk1 Þ þ vjlhk1
ð90Þ
¼ f T ðxjlhk1 Þ þ vjlhk1 .
ð91Þ
From Eqs. (86) to (87), a partition of N samples into N=2L blocks is used where L is the level number of the multires decomposition. Since it is impossible to analytically derive f T ðÞ in Eq. (91), we take another approach. ~ jlh , we can find the For the thresholded o k1 j corresponding x~ lhk1 . Therefore, we can just propagate the elements in x~ jk1 ¼ TT x~ jlhk1 which corre~ jk1 ¼ TT o ~ jlh . spond to the distinct elements in o k1
1 ~ ihk1 . ¼ pðzk jfxhjm gikjk1 Þo C
Nl X
ð85Þ
N=2L
sgnðfxhjm gikjk1 Þ
hj XX X
(84)
and the propagated density is described by
N hj
þ
vik Nð0; Qk Þ
m
~ ihk fxhjm gikjk1 o ð83Þ
~ jk1 , we can For the repeated elements of o propagate one representative element in x~ jk1 .The number of particle propagations is controlled by the thresholding operation, which can significantly reduce the computation in propagation while maintaining performance. The details of finding the corresponding elements in x~ jlhk1 given ~ jlh are described as follows. For an example o k1 of three levels (L ¼ 3), we can explicitly derive
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
the form of 2 3 h 6 3 6 h 6 6 2 6 h 6 6 6 0 6 T¼6 6 h 6 6 6 0 6 6 6 0 4
~ jlh ð2Þ ¼ 0()sum of 1st 4 and sum of 2nd 4 o k1 ~ jk1 are the same; elements of o j ~ lh ð3Þ ¼ 0()sum of 1st 2 and sum of 2nd 2 o k1 ~ jk1 are the same; elements of o j ~ lh ð4Þ ¼ 0()sum of 3rd and sum of 4th 2 o k1 ~ jk1 are the same; elements of o j ~ lh ð5Þ ¼ 0()1st and 2nd elements of o ~ jk1 o k1 are the same; ~ jlh ð6Þ ¼ 0()3rd and 4th elements of o ~ jk1 o k1 are the same; ~ jlh ð7Þ ¼ 0()5th and 6th elements of o ~ jk1 o k1 are the same; ~ jlh ð8Þ ¼ 0()7th and 8th elements of o ~ jk1 o k1 are the same.
T as h3
h3
h3
h3
h3
h3
h3
h3
h3
h3
h3
h3
h2
h2
h2
0
0 2
h2
0
0
0
h
h
0
0
0
0
0
0
h
h
0
0
0
0
0
0
h
h
0
h
7 h3 7 7 7 07 7 7 h2 7 7 7, 07 7 7 07 7 7 07 5
0 2
3
h3
0 0 0 0 0 0 h h pffiffiffi 2 h¼ . ð92Þ 2 Here is a list of correspondences between zero ~ jlh ~ jk1 : elements in o and the same elements in o k1
1397
A combination of the above correspondences could also result in a combination of the same ~ jlh ð2Þ ¼ 0 and element numbers. For instance, fo k1 j ~ lh ð7Þ ¼ 0g(){sum of 1st 4 and sum of 2nd 4 o k1 ~ jk1 are the same, and 5th and 6th elements of o
~ jlh ð1Þ ¼ 0()sum of 8 elements of o ~ jk1 is o k1 zeros;
original p (x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x unires propagated p (x1|0)
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
x unires updated p (x1|1)
p (x1|1)
0.03 0.02 0.01 0 -5
0
5
10
15
20
25
x Fig. 11. Traditional (unires) particle filtering with 1000 particles: (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1398
~ jk1 rm are the same} or fo ~ jlh ð2Þ ¼ elements of o k1
where wk Nð0; Rk Þ, we want to update the propagated density
~ jlh ð5Þ ¼ 0 and ojlh ð6Þ ¼ 0g(){sum of 1st 4 0; o k1 k1 ~ jk1 are the same, 1st and sum of 2nd 4 elements of o j ~ k1 are the same and 3rd and and 2nd elements of o ~ jk1 are the same}. 4th elements of o The propagation of the thresholded quantity is then given by x~ jkjk1 ¼ f ðTT x~ jlhk1 Þ þ vjk1
pðxk jz1:z1 Þ ¼
(93)
~ jk1 . dT ðxjk x~ jkjk1 Þo
l i ¼ Nðzk gðx~ ikjk1 Þ; Rk Þ;
~ i ¼ 1; . . . ; N,
(94)
Explicit multires update: Given the measurement model
pðxk jz1:z Þ ¼
zk ¼ gðxk Þ þ wk ,
where C ¼
(95)
N~ X 1 i i ~ k1 l dðxk x~ ikjk1 Þ, o C i¼1
PN~
~ ik1 i¼1 o
li .
original unires p (x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x
b
multires propagated p (x1|0), thresholded
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
15
20
x
c
multires updated p (x1|1), thresholded
p (x1|1)
0.04 0.02 0 -5
0
(97)
where N~ is the total number of distinct particles, which could be much smaller than N depending on the threshold. The density update can be carried out by
j¼1
a
(96)
with the explicit multires method. The likelihood functions are
N=2L
X
~ jk1 dT ðxjk x~ jkjk1 Þo
j¼1
which in turn defines the density propagation as pðxk jz1:k1 Þ ¼
L N=2 X
5
10 x
Fig. 12. Multires particle filtering with 562 particles (T s ¼ 0:0002): (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
(98)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1399
Example 4. The nonlinear system and measurement models are given by
without the multiresolution approach, where
xk1 5xk1 þ xk ¼ þ vk ; 2 ð1 þ x2k1 Þ
x^ 0j0 ¼
vk Nð0; 1Þ
(99)
zk ¼ 10 tan1
k
10
þ wk ;
xi0 oi0
i¼1
and x
N X
wk Nð0; 1Þ
with xi0 px ðx0 Þ. Fig. 11 shows the propagation and update of px ðx0 Þ using the traditional particle filter approach (we call it the unires approach) with the number of particles being 1000. It can be seen that the particle cloud peak moves from 10 to 6 with a wider spread in (b) due to propagation and becomes more concentrated in (c) due to update. Figs. 12 and 14 show the propagation and update of px ðx0 Þ with the newly developed multires particle approach (only the explicit approach is used due to the paper length limit) where the threshold (T s ) is 0.0002 in Fig. 12, 0.001 in Fig. 13 and 0.002 in Fig. 14. The propagation-update cycle from t ¼ 0 to 1 can be applied to other time instants in any
(100)
with the initial nonlinear (bi-modal) state density being px ðx0 Þ ¼ 0:7 Nð10; 22 Þ þ 0:3 Nð20; 42 Þ.
(101)
In the following, we will show how px ðx0 Þ is propagated from t ¼ 0 to 1 using particles with and without the multiresolution approach. Then the propagated density will be updated by z1 ¼ 15 tan1 ðx^ 0j0 =10Þ þ w1 using particles with and
a
original unires p (x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x
b
multires propagated p (x1|0), thresholded
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
15
20
x
c
multires updated p(x1|1), thresholded
p (x1|1)
0.04 0.02 0 -5
0
5
10 x
Fig. 13. Multires particle filtering with 234 particles ðT s ¼ 0:001Þ: (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1400
a
original unires p(x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x
b
multires propagated p(x1|0), thresholded
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
x
c
multires updated p (x1|1), thresholded
p (x1|1)
0.04 0.02 0 -5
0
5
10
15
20
x Fig. 14. Multires particle filtering with 165 particles ðT s ¼ 0:002Þ: (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
Table 1 A comparison of computational complexity and estimation accuracy (500 runs) Approaches
Unires particle filter Multires particle filter ðT s ¼ 0:0002Þ Multires particle filter ðT s ¼ 0:001Þ Multires particle filter ðT s ¼ 0:002Þ
RMS estimation errors
Number of particles
0.0450 0.0451
1000 562
0.0443
234
0.0452
165
filtering application; therefore, a complete multires particle filtering algorithm is derived. A numerical comparison of estimation accuracy and computational complexity of both unires and multires approaches is given in Table 1. It can be clearly seen
that the proposed multires approach performs comparably (or even slightly better sometimes) to the unires approach, but is significantly more efficient! 4. Conclusions We have developed a framework for computationally efficient particle filtering using the spatial domain multiresolutional approach. Both implicit and explicit implementation methods for a multires particle filter are presented. Computational efficiency of a particle filter is achieved by drastically reducing the needed number of particles while maintaining a comparable estimation accuracy. References [1] N.J. Gordon, D.J. Slamond, A.F.M. Smith, Novel approach to nonlinear/non-Gaussian Bayesian state estimation, IEE Proc.-F 140 (2) (April 1993) 107–113.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401 [2] R. Van der Merwe, N. de Freitas, A. Doucet, E.A. Wan, The unscented particle filter, Technical Report CUED/F-INFENG/TR380, Cambridge University Engineering Department, August 2000. [3] A. Doucet, N. de Freitas, N. Gordon, Sequential Monte Carlo Methods in Practice, Springer, New York, 2001. [4] N. Cui, L. Hong, J.R. Layne, A comparison of nonlinear filtering approaches with application to ground target tracking, Signal Processing 85 (8) (August 2005) 1469–1492. [5] A. Farina, B. Ristic, Tracking a ballistic target: comparison of several nonlinear filters, IEEE Trans. Aerospace Electronic Systems 38 (3) (July 2002) 854–867. [6] A. Doucet, N. Gordon, Efficient particle filters for tracking maneuvering targets in clutter, in: IEE Colloquium on Target Tracking: Algorithms and Applications (Ref. No. 1999/090, 1999/215), 1999, pp. 4/1–4/5.
1401
[7] D. Fox, KLD-sampling: adaptive particle filter, Adv. Neural Inf. Process. Syst. (NIPS) 14 (2001) 713–720. [8] L. Hong, Multiresolutional distributed filtering, IEEE Trans. Automat. Control 39 (4) (1994) 853–856. [9] L. Hong, Multiresolutional multiple-model target tracking, IEEE Trans. Aerospace Electronic Systems 30 (2) (April 1994) 518–524. [10] L. Hong, Multirate interacting multiple model filtering for target tracking using multirate models, IEEE Trans. Automat. Control 44 (7) (1999) 1326–1340. [11] L. Hong, N. Cui, M. Bakich, J.R. Layne, Multirate interacting multiple model particle filter with application to terrain-based ground target tracking, IEE Proc: Control Theory Appl. 153 (6) (2006) 721–731. [12] I. Daubechies, Orthonormal bases of compact supported wavelets, Commun. Pure Appl. Math. XLI (1988) 909–996.
Signal Processing 87 (2007) 1384–1401 www.elsevier.com/locate/sigpro
A spatial-domain multiresolutional particle filter with thresholded wavelets Lang Honga,, Devert Wickerb a
Department of Electrical Engineering, Wright State University, Dayton, OH 45435, USA Target Recognition Branch, Air Force Research Laboratory, SNAT, WPAFB, OH 45433, USA
b
Received 10 July 2006; received in revised form 19 July 2006; accepted 5 December 2006 Available online 10 January 2007
Abstract Particle filters have been proven to be very effective for nonlinear/non-Gaussian filtering. However, the most notorious disadvantage of a particle filter is its formidable computational complexity, since hundreds (even thousands) of particles are usually needed to achieve a required approximation accuracy. It has also been proven that one of the techniques of truly solving a computational problem is multiresolutional processing, both in temporal and spatial domains. Therefore, in this paper we propose a multiresolutional particle filter in the spatial domain using thresholded wavelets to reduce significantly the number of particles, meanwhile maintaining the full strength of a particle filter. r 2007 Elsevier B.V. All rights reserved. Keywords: Particle filters; Nonlinear/non-Gaussian estimation; Computational efficiency; Multiresolutional processing techniques; Wavelets
1. Introduction Particle filters, falling into the category of sequential Monte Carlo methods, are effective in nonlinear and non-Gaussian filtering. A particle filter estimates a posterior density function pðxk jZ k Þ by an empirical approach ^ k jZ k Þ ¼ pðx
N 1X dx ðxkðiÞ Þ, N i¼1
(1)
where xkðiÞ are iid (independently and identically distributed) random samples of pðxk jZ k Þ and Zk represents the accumulative measurements. The Corresponding author.
E-mail addresses: [email protected] (L. Hong), [email protected] (D. Wicker).
associated integration can also be achieved empirically as Z N 1X f ðxðiÞ Þ. (2) Iðf t Þ ¼ f t ðxk Þpðxk jZ k Þ dxk N i¼1 t k Although operations in Eqs. (1) and (2) are straightforward and can be applied to any nonlinear and non-Gaussian density, there is a problem: sampling directly from a posterior function is not always possible. Researchers have found a solution: instead of sampling the posterior function, samples can be drawn from a proposal function whose support includes that of the posterior function. This method is called importance sampling (IS) or sequential importance sampling (SIS) when the weight updates are performed in a recursive manner. However, IS/SIS is not problem free. As time elapses, the
0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2006.12.005
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
distribution of importance weights becomes more and more skewed. After a few time steps, only a few particles have non-zero weights. Another fix has also been proposed: a selection step called importance resampling is introduced to eliminate the particles with low weights and multiply the particles having high weights. More recently, many improvements have been made and the technology of particle filters is becoming mature, which includes bootstrap particle filters, extended-Kalman particle filters, unscented particle filters, etc. [1,2]. Ref. [3] presents an excellent collection of articles on particle filters. When applied to target tracking, both references [4,5] give a good comparison of different nonlinear filtering approaches, including particle filters. Although particle filters are effective for nonlinear/non-Gaussian filtering, the most notorious disadvantage of a particle filter is its formidable computational complexity, since hundreds (even thousands) of particles are usually needed to achieve a decent density approximation accuracy. Little effort has been attempted with limited success on improving the computational efficiency of a particle filter [4,6,7]. It is well known that there are two fundamental approaches in achieving computational efficiency in engineering applications: parallel processing and multiresolutional processing. Parallel processing divides a computational task into many subtasks which can then be executed simultaneously with multiple processors, while multiresolutional processing prioritizes the computational resources according to the need to maintain the ‘‘details’’ of signals. Many engineering applications cannot be implemented in parallel due to the nature of a single computational thread. On the other hand, multiresolutional processing has a broader application range. This paper presents a computationally efficient particle filter using a wavelet-based multiresolutional estimation technique developed over the last decade by the first author [8–10]. In general, multiresolutional processing can be applied to a particle filter in either temporal or spatial domains. A multiresolutional particle filter is called a multirate particle filter if it is derived by applying multiresolutional processing in a temporal domain. On the other hand, applying multiresolutional processing in a spatial domain in designing a particle filter results in a multires particle filter. In our previous work, we have demonstrated significant computational savings by using a multirate multiple model particle filter [11]. This paper is focused on computational savings with a multires
1385
particle filter. A particle filter approximates a density function by a set of particles with weights and a multires particle filter further decomposes the weights into lowpassed and highpassed components. By thresholding the highpassed particle weights, a much smaller number of particles can represent the density function with a comparable approximation accuracy. A density can be approximated by a set of densities with different resolutions (smoothness), where a smoother density requires a smaller number of particles. As a result, the total number of particles could be much smaller than a regular particle filter (we call this a unires particle filter), while maintaining a comparable approximation accuracy. The rest of the paper is organized as follows. Section 2 formulates a multires particle filter and Section 3 presents detailed implementations of a multires particle filter with simulation examples. Section 4 concludes the paper. 2. Formulation of a multires particle filter For a nonlinear filtering problem, the system and measurement models are given by xk ¼ f ðxk1 Þ þ vk1 ,
ð3Þ
zk ¼ hðxk Þ þ wk ,
ð4Þ
where f ðÞ and hðÞ are system dynamic and measurement functions. Variables vk1 and wk are governed by the following known densities, pv and pw : pðxk jxk1 Þ ¼ pv ðxk f ðxk1 ÞÞ, pðzk jxk Þ ¼ pw ðzk hðxk ÞÞ.
ð5Þ ð6Þ
The initial state density, px ðx0 Þ, is also assumed known which can be decomposed into a multires structure with j 1 levels X XX px ðx0 Þ ¼ aj1 ;k fj1 ;k ðx0 Þ þ bj;k cj;k ðx0 Þ (7) k
jXj 1
k
with scaling and wavelet coefficients given by [12] aj1 ;k ¼ hpx ðx0 Þ; fj 1 ;k ðx0 Þi,
ð8Þ
bj;k ¼ hpx ðx0 Þ; cj;k ðx0 Þi,
ð9Þ
where ffj1 ;k ¼ 2j 1 =2 fð2j1 x kÞ; k 2 Zg, fcj;k ¼ 2j=2 cð2j x kÞ; k 2 Z; jXj 1 g, and fðxÞ and cðxÞ are the scaling and mother wavelet functions, respectively. In the form of
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1386
particle representation, we have px ðx0 Þ ¼
N0 X
oix0 dðx0 xi0 Þ,
(10)
i¼1
where particle xi0 is sampled from px ðx0 Þ, i.e., xi0 px ðx0 Þ and oix0 is the corresponding weight. The inner products defined by Eqs. (8) and (9) that require nonlinear integration become a^ j1 ;k ¼
N0 X
ol 1 ¼
N 0 =2 X
i
(18)
i
(19)
oil 1 dðx0 x01 Þ
i¼1
and oix0 fj1 ;k ðxi0 Þdðx0
xi0 Þ
(11)
i¼1
oh1 ¼
and b^ j;k ¼
where only level 1 decomposition is presented. Not considering the edge effect, the process of convolution and downsampling results in N20 non-zero coefficients which can be described in a summation form
N 0 =2 X
oil 1 dðx0 x01 Þ,
i¼1 N0 X
i
oix0 cj;k ðxi0 Þdðx0
xi0 Þ.
(12)
i¼1
In this paper, we will use Harr wavelets as a vehicle for demonstration.1
where x01 is the result of lowpassed and downsampled xi0 . Level 2 decomposition can be derived from ol 1 by ol 2 ½i ¼
þ1 X
ol 1 ½i h½2i k
(20)
ol 1 ½i g½2i k,
(21)
k¼1
2.1. Harr wavelets
and
The scaling function is a simple rectangular function 1 if 0ptp1; fðtÞ ¼ (13) 0 otherwise with only two non-zero coefficients hð0Þ ¼ hð1Þ ¼ p1ffiffi2 and the wavelet function is 8 if 0pto0:5; > <1 (14) cðtÞ ¼ 1 if 0:5ptp1; > : 0 otherwise with two non-zero coefficients gð0Þ ¼ p1ffiffi2 and gð1Þ ¼ p1ffiffi2. A filterbank is used for the discrete wavelet transform (DWT). Rewrite Eq. (10) as a signal sequence px ðx0 Þ ¼ fox0 ½ig;
i ¼ 1; . . . ; N 0 .
(15)
Eqs. (11) and (12) can be expressed as a discrete convolution with a downsampling process ol 1 ½i ¼ a^ 1 ½i ¼
þ1 X
ox0 ½i h½2i k
(16)
oh2 ½i ¼
þ1 X k¼1
where only N40 non-zero coefficients are obtained. Other levels of decomposition can be obtained similarly. The major difference from the traditional wavelet transform in Eqs. (11) and (12) is that xi0 ; i ¼ 1; . . . ; N are random samples, not a predefined and evenly spaced sequence. To overcome the randomness, we sample all N 0 particles first; then filtering is applied, but the samples are still unevenly spaced. Because of the downsampling operation, the numbers of transformed particles ol j1 and ohj are drastically reduced with the number of levels. With the scaling and wavelet coefficients estimated by Eqs. (16) and (17), the density (weights) reconstruction can be performed by ox0 ½i ¼
þ1 X
ol 1 ½kh0 ½2k i þ oh1 ½kg0 ½2k i.
k¼1
k¼1
(22)
and oh1 ½i ¼ b^ 1 ½i ¼
þ1 X
ox0 ½i g½2i k,
(17)
k¼1 1
The multires technology developed here can be applied to any wavelet.
A simple thresholding can be carried out on the wavelet coefficients to reduce the number of coefficients. Let ~ h1 ½i ¼ oh1 ½i if oh1 ½iXT s ; o
~ h1 ½i ¼ 0. otherwise o
(23)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
A thresholded density estimation can then be achieved as þ1 X
ox0 ½i ¼
~ h1 ½kg0 ½2k i. ol 1 ½kh0 ½2k i þ o
k¼1
(24) T s is an application dependent parameter which needs to be adjusted for each application depending on available computational capability. The larger T s is, the fewer number of surviving coefficients but larger approximation errors. Example 1. Assume density px ðx0 Þ is given by the following bi-modal function: px ðx0 Þ ¼ 0:6Nð10; 22 Þ þ 0:4Nð20; 42 Þ,
(25)
where NðÞ denotes a normal distribution function. px ðx0 Þ is approximated by Eq. (10) whose particle weights are plotted in Fig. 1 with the number of particles being 1000. A three-level multires decomposition for the particle weights is performed and the results are shown in Fig. 2. The original particle weights are now represented by multiresolutional weights in Fig. 2, where (a) and (b) are level-1 lowpass-filtered and highpass-filtered weights, (c) and (d) are level-2 lowpass-filtered and highpass-filtered weights, and (e) and (f) are level-3 lowpass-filtered and highpass-filtered weights. In the reconstruction process, weights from (b), (d)–(f) completely reconstruct the weights in Fig. 1. When a threshold is applied to the highpassed particle weights, the numbers of particles can be drastically reduced at each level without significantly degrading reconstruction accuracy. Figs. 3 and 4 show the thresholding results, where Fig. 3 shows that only a small number of highpass-filtered particle weights
are used in the reconstruction with the threshold being 0.001. The results for threshold 0.002 are given in Figs. 5 and 6. Thus far, the decomposition is performed for the sequence where the independent variable is the sequence index i. If the independent variable is a continuous variable, such as x, the decomposition is handled differently. Assume density px can be approximated by px ¼
N X
oi dðx xi Þ,
(26)
i¼1
where xi is sampled from px , i.e., xi px , and P N i¼1 oi ¼ 1. Density px is a continuous function of x, not i. To better handle Eq. (26) in multiresolution decomposition, we can rewrite it in a vector form px ¼ dT o ,
(27)
where d ¼ ½dðx x1 Þ; . . . ; dðx xN ÞT , o ¼ ½o1 ; . . . ; oN T .
ð28Þ
Let o ¼ TT ½ol j ; ohj ; . . . ; ohj 1
1
1 þL
px ¼ dT TT olh ¼ ðT d ÞT olh ¼ dTlh olh ,
0.008 0.006 0.004 0.002 0 0
5
10
15
20
25
30
35
x Fig. 1. The weights of a particle approximation of the density defined in Eq. (25).
(29)
(30)
where 1
0.01
T ¼ TT olh
be the multiresolution decomposition of sample weights, where T is an orthogonal wavelet decomposition transform matrix, such that T1 ¼ TT . There are L þ 1 levels in the decomposition. Substituting Eq. (29) into Eq. (27) results in
dlh ¼ T d ¼ ½dl j ; dhj ; . . . ; dhj
0.012
particle weights
1387
1
1 þL
T .
(31)
The multiresolution decomposition in Eq. (30) is applied to both weights and the independent variable.2 There are two methods of implementing Eq. (31): implicit and explicit. The implicit method embeds the wavelet transformation into a complicated variable structure and no inverse transform is needed when performing density reconstruction. On the other hand, the explicit method uses a simple variable structure, but an inverse transform is needed in density reconstruction. 2
To be precise, multiresolution decomposition is applied to the d function of the independent variable.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1388
b
a
x 10-3 4
0.008
level 1 highassed
level 1 lowpassed
0.01
0.006 0.004 0.002 0
2 0 -2 -4
0
5
10
15
20
25
30
35
0
5
10
15
x
c
d
x 10-3
level 2 highpassed
level 2 lowpassed
25
30
35
20
25
30
35
20
25
30
35
x 10-3 4
8 6 4 2 0
2 0 -2 -4
0
5
10
15
20
25
30
35
0
5
10
15
x
e
x
f
x 10-3
1 level 3 highpassed
6 level 3 lowpassed
20 x
4
2
0
x 10-3
0.5 0 -0.5 -1 -1.5
0
5
10
15
20
25
30
35
0
5
10
x
15 x
Fig. 2. A three-level multires decomposition of the particle weights.
Implicit method: To make an equivalent variable for each dl and dh in Eq. (31), we introduce the following variable structure:
and sgn is a sign function ( 1; x 2 first half of fxhj g; sgnðfxhj gÞ ¼ 1; x 2 second half of fxhj g:
dl j ¼ hj 1 dðx fxl gÞ,
(34) With this implicit method, the density reconstruction can be performed as follows: Nl X px ¼ oil dðx fxl gi Þ
1
pffiffiffi 2 j , ð32Þ dhj ¼ h dðx fxhj gÞsgnðfxhj gÞ; jXj 1 ; h ¼ 2 where fxl g and fxhj g are collections of sample points
i N
fxl g ¼ fsamples associated with ol g fxhj g ¼ fsamples associated with ohj g
and
þ ð33Þ
hj XX
jXj 1
i
oihj dðx fxhj gi Þsgnðfxhj gi Þ.
ð35Þ
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
a
b
x 10-3 4
level 1 highassed
level 1 lowpassed
0.01 0.008 0.006 0.004 0.002 0
2 0 -2 -4 -6
0
5
10
15
20
25
30
35
0
5
10
15
x
c
d level 2 highpassed
level 2 lowpassed
6 4 2 0
30
35
20
25
30
35
20
25
30
35
x 10-3 2 1 0 -1 -2
0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
x 10-3 6
x 10-3 2
level 3 highpassed
level 3 lowpassed
25
3
8
4
2
0
20 x
x 10-3
e
1389
0
5
10
15
20
25
30
35
1 0 -1 -2
0
5
10
15
x
x
Fig. 3. A three-level multires decomposition of the particle weights with thresholding where the threshold is 0.001.
Explicit method: It can be seen in Eq. (35) that no explicit inverse transform is needed in density reconstruction, but the variable structure in Eq. (32) may pose a difficulty in implementation. Other than a complicated variable structure, we can introduce transformed variables as 3 2 2 13 xl x 6 x 7 6 x2 7 6 h1 7 6 7 7 6 7 (36) xslh ¼ 6 .. 7 ¼ T6 6 .. 7. 6 . 7 4 5 . 5 4 xhL2L xN
Then by using dlh ¼ T d ¼ Tdðx TT xslh Þ,
(37)
we have px ¼ ðTdðx TT xslh ÞÞT olh T
¼ ðTdðT ðxlh
xslh ÞÞÞT olh
ð38Þ ð39Þ
which contains explicit transformation T. Eq. (38) can be rewritten as px ¼ pT ðdlh ; olh Þ, where pT is a transformation mapping function.
(40)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1390
b
0.01 0.008
original level 2 lowpassed
original level 2 lowpassed
a
0.006 0.004 0.002
0.01 0.008 0.006 0.004 0.002
0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
d
0.01 0.008
original level 1 lowpassed
original level 1 lowpassed
c
0.006 0.004 0.002
30
35
20
25
30
35
20
25
30
35
0.008 0.006 0.004 0.002 0
0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
0.01 0.008
0.01
original level 0 lowpassed
original level 0 lowpassed
25
0.01
0
e
20 x
0.006 0.004 0.002
0.008 0.006 0.004 0.002
0
0 0
5
10
15
20
25
30
35
0
5
x
10
15 x
Fig. 4. A density reconstruction with thresholding where the threshold is 0.001.
If we name each level’s density component3 in Eq. (30) by pjx1l ¼ dTl ol
(41)
for the lowpassed density component at level j 1 for the implicit method or pjx1l ¼ pT ðdl ; ol Þ
(42)
for the explicit method, and pjxh ¼
X
dThjk ohjk ;
j ¼ j1; . . . ; L
(43)
k 3 We call them density components, since they are not regular densities anymore, due to the fact that the integration of each component is not one.
for the highpassed density component at level j using the implicit method, or X pT ðdhjk ; ohjk Þ; j ¼ j 1 ; . . . ; L (44) pjxh ¼ k
using the explicit method, we could rewrite Eq. (30) as X px ¼ pjx1l þ pjxh . (45) jXj 1
With simple thresholding, we can significantly reduce the number of particles without sacrificing the performance. Thresholding can be applied to both the implicit and explicit methods. For the implicit method, thresholding is applied to the density components 8 T < fdhjk ohjk gi ; jfdThjk ohjk gi j4T s ; T i ~ hjk g ¼ fdhjk o (46) : 0; jfdThjk ohjk gi jpT s
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
a
b
x 10-3 4
level 1 highpassed
level 1 lowpassed
0.015 0.01 0.005 0
2 0 -2 -4
0
5
10
15
20
25
30
35
0
5
10
15
x
c
d
30
35
20
25
30
35
20
25
30
35
5 level 2 highpassed
level 2 lowpassed
25
x 10-3
8 6 4 2
4 3 2 1
0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
x 10-3
1 level 3 highpassed
8 level 3 lowpassed
20 x
x 10-3
e
1391
6 4 2
0.5 0
-0.5
0
-1 0
5
10
15
20
25
30
35
0
5
10
15
x
x
Fig. 5. A three-level multires decomposition of the particle weights with thresholding where the threshold is 0.002.
and for the explicit method, thresholding is directly applied to the wavelet coefficients ( ohjk ; johjk j4T s ; ~ hjk ¼ o (47) 0; johjk jpT s : Thus, we have X ~ hjk ; dThjk o pjxh ¼
Example 2. Let the wavelet transform be 3 2 1=2 1=2 1=2 1=2 7 6 1=2 1=2 7 6 1=2 1=2 pffiffiffi pffiffiffi 7, T¼6 6 2=2 2=2 0 0 7 4 pffiffiffi pffiffiffi 5 0 0 2=2 2=2
(50)
and a 4-point density approximation be j ¼ j1 ; . . . ; L
(48)
k
for the implicit method, and X ~ hjk Þ; j ¼ j 1 ; . . . ; L pjxh ¼ pT ðdhjk ; o k
for the explicit method.
½x1 x2 x2 x4 ¼ ½2 10 4 6T , ½o1 o2 o3 o4 ¼ ½16
(49)
1 5 1T 4 12 6
ð51Þ
and 5 px ¼ 16 dðx 2Þ þ 12 dðx 4Þ þ 16 dðx 6Þ þ 14 dðx 10Þ.
(52)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1392
b
0.01
0.01
original level 2 lowpassed
original level 2 lowpassed
a
0.008 0.006 0.004
0.008 0.006 0.004 0.002
0.002 0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
d
0.01
original level 1 lowpassed
original level 1 lowpassed
c
0.008 0.006 0.004
30
35
20
25
30
35
20
25
30
35
0.008 0.006 0.004 0.002
0.002
0 0
5
10
15
20
25
30
35
0
5
10
15
x
x
f
0.01
0.01
original level 0 lowpassed
original level 0 lowpassed
25
0.01
0
e
20 x
0.008 0.006 0.004
0.008 0.006 0.004 0.002
0.002
0
0 0
5
10
15
20
25
30
35
0
5
10
15
x
x
Fig. 6. A density reconstruction with thresholding where the threshold is 0.002.
Implicit method: From Eqs. (29) and (31), we have 2
dl
3
7 6 6 dh 1 7 7 6 7 ¼ Td 6 6 dh21 7 5 4 dh22 2 3 1=2½dðx 2Þ þ dðx 4Þ þ dðx 6Þ 6 7 6 7 þdðx 10Þ 6 7 6 7 6 1=2½dðx 2Þ dðx 4Þ þ dðx 6Þ 7 7 6 ¼6 7 7 6 þdðx 10Þ 7 6 pffiffiffi 7 6 7 6 2 =2½dðx 2Þ þ dðx 4Þ 5 4 pffiffiffi 2=2½dðx 6Þ þ dðx 10Þ ð53Þ
and 2 1 3 3 2 ol 2 6 17 6o 7 6 12 7 6 h1 7 6 pffiffi 7 7 6 7. 6 oh 7 ¼ T o ¼ 6 6 82 7 4 21 5 4 pffiffi 5 2 oh22
(54)
24
Then px ¼ dl ol þ dh1 oh1 þ dh21 oh21 þ dh22 oh22 ð55Þ 1 ¼ 4 ½dðx 2Þ þ dðx 4Þ þ dðx 6Þ þ dðx 10Þ 1 24 ½dðx 2Þ dðx 4Þ þ dðx 6Þ þ dðx 10Þ þ 18 ½dðx 2Þ 1 þ dðx 4Þ þ 24 ½dðx 6Þ þ dðx 10Þ. ð56Þ
For this example, fxl g ¼ f2; 4; 6; 10g, fxh1 g ¼ f2; 4; 6; 10g, fxh21 g ¼ f2; 4g and fxh22 g ¼ f6; 10g. If we set
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1393 j
1 the threshold as 24 , applying thresholding to Eq. (56) results in
Figs. 7 and 8 show the reconstruction of px from px1l and pjxh with and without thresholding using the implicit method.
px dl ol þ dh1 oh1 þ dh21 oh21 þ dh22 oh22
Explicit method: Using Eq. (36), we can find the following transformed variables: pffiffiffi pffiffiffi xlh ¼ ½11 5 2 2 2T . (59)
¼
1 4 ½dðx 2Þ þ dðx 4Þ þ dðx þ 18 ½dðx 2Þ þ dðx 4Þ.
ð57Þ
6Þ þ dðx 10Þ ð58Þ
5 5 12 12
px
Threshold = 0 4 multires particles reconstruct original 4 particles
1 1 6 6
1 1 4 4
1 1 6 6
x
px
h2
2 1 24
6 4
- 1 8
6
8
10
4 particles
1 8 2
4
8 - 1 24 2 particles
Direct sum; no inverse transformation is needed; but complicated variable structure is needed
+
x
10
px
px
l
h1
1 4 1 24 2
6 4
8
10
x
x
- 1 24
4
2
6
1 particle
8
10
1 particle j
Fig. 7. Multires density reconstruction from px1l and pjxh using the implicit method without thresholding.
px Threshold = 1 24 2 multires particles approximate original 4 particles
1 6 1 8
5 12 3 8
x
px
2
h2
2
- 1 8
8
10
+ Direct sum; no inverse transformation is needed; but complicated variable structure is needed
x 6 1 particle
8 6 4 particles
4
1 8
4
1 1 4 4
1 1 4 6
10
px
px
l
h1
1 4 x 2
4
6 0 particle
8
x
10
2 j
4
6 1 particle
8
10
Fig. 8. Multires density reconstruction from px1l and pjxh using the implicit method with thresholding.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1394
The reconstruction of the density can be performed separately by variable and weight reconstructions, which are
a density in a multires decomposition form X pjxh , pðxk1 jz1:k1 Þ ¼ pjx1l þ
x ¼ ½2 4 6 10T ¼ TT xlh
where
(60)
and o¼
1
1 T
5 1 6 12 6 4
pjx1l
T
¼ T olh .
(61)
The density reconstruction is then px ¼
1 6 dðx
2Þ þ
k1
5 12 dðx
4Þ þ
1 6 dðx
6Þ þ
1 4 dðx
10Þ
whichpisffiffi illustrated in Fig. 9. If the threshold is set as T s ¼ 242, then we have pffiffi T 2 1 olh ¼ ½12 12 (63) 8 0 ,
¼ dTl k1 ol k1
and pjxh
(62)
¼ k1
X m
dThjm
(66)
k1
~ hjm o
k1
;
j ¼ j 1 ; . . . ; L.
(67)
In the following, we will derive a complete update cycle for multires densities. 3.1. The implicit method Multires propagation:
which results in o ¼ ½16
k1
(65)
k1
jXj 1
5 5 5 T 12 24 24 .
(64)
Figs. 9 and 10 demonstrate the results of density reconstruction using the explicit method, with and without thresholding. 3. Multires bootstrap particle filter algorithm
pðxk jz1:k1 Þ Z ¼ pðxk jxk1 Þpðxk1 jz1:k1 Þ dxk1 Z ¼ pðxk jxk1 Þpjx1l dxk1 Z X þ pðxk jxk1 Þ pjxh dxk1 ¼
Nl Z X
pðxk jfxl gik1 Þoil k1 dðxk1 fxl gik1 Þ dxk1
i¼1
5 5 12 12
px
1 1 4 4
1 1 6 6
1 1 6 6
Threshold = 0 4 multires particles reconstruct original 4 particles
x
ωh2
24
8
6
10
4 particles 2 24
TT
xh2
2 22 2 4
6
8
Need an inverse transform
10
2 particles 1 2
ωl
ωh1
xh1
5 2
4 6 8 1 1 particle 12
10
ð69Þ
jXj 1
This section presents a major development of a spatial domain multiresolution (multires) particle filter for computationally efficient processing. Given
2 8
ð68Þ
xl 2
4
6
8
10 11
1 particle
Fig. 9. Multires density reconstruction from olh and xlh using the explicit method without thresholding.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
5 5 12 12
px 2 Threshold = 24 3 multires particles reconstruct original 4 particles
1395
1 6
1 4
1 5 6 24
1 6
5 24 x
ωh2 2
4
8
6
10
4 particles
2 8
TT
xh2 4
2 2
6 8 1 particle
Need an inverse transform
10
1 2
ωl
ωh1
xh1
5 2
4 6 8 - 1 12 1 particle
10
xl 2
4
6 1 particle
8
10 11
Fig. 10. Multires density reconstruction from olh and xlh using the explicit method with thresholding.
hj XZ XX
N
þ
jXj 1 i¼1
m
The propagated state estimate can be derived as pðxk jfxhjm gik1 Þoihk1
dðxk1 fxhjm gik1 Þsgnðfxhjm gik1 Þ dxk1 ¼
Nl X
x^ kjk1 ¼ ð70Þ
Nl X
pðxk jfxl gik1 Þoil k1
þ
jXj 1 i¼1
m
~ ihk1 pðxk jfxhjm gik1 Þo
sgnðfxhjm gik1 Þ ¼
Nl X
oil k1 dðxk
ð71Þ
2=3
jXj 1 i¼1
m
vk Nð0; s2v Þ, 2=3
~ ihk1 o
sgnðfxhjm gikjk1 Þ,
dðxk
ð73Þ
(74)
where sv ¼ 1. The propagation density is pðxk jxk1 Þ ¼ Nðxk xk1 ; s2v Þ.
N
þ
~ ihk1 fxhjm gikjk1 o
Example 3. Continuing with Example 2. Assume the system model is xkþ1 ¼ xk þ vk ;
fxl gikjk1 Þ
i¼1 hj X XX
m
sgnðfxhjm gikjk1 Þ.
N
þ
N hj jX 1 þL X X j¼j 1 i¼1
i¼1 hj X XX
oil k1 fxl gikjk1
i¼1
fxhjm gikjk1 Þ ð72Þ
where fxl gikjk1 pðxk jfxl gik1 Þ and fxhjm gikjk1 pðxk j fxhjm gik1 Þ. Since both fxl gikjk1 and fxhjm gikjk1 are associated with a same set of particles at the original resolution level, only on one propagated set (either fxl gikjk1 or fxhjm gikjk1 ) is obtained from sampling and others can be derived from the sample set.
(75)
Take fxl gk1 ¼ f2; 4; 6; 10g for propagation. By sampling pðxk jfxl gik1 Þ, we have fxl gkjk1 ¼ f1:5; 2:8; 3:5; 5:4g. Due to the relationship defined by T in Eq. (50), we can derive other propagated sets fxh1 gkjk1 ¼ f1:5; 2:8; 3:5; 5:4g;
and
fxh22 gkjk1 ¼ f3:5; 5:4g.
fxh21 gkjk1 ¼ f1:5; 2:8g
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1396
Multires update: The update of propagated multires density is given by
There will be a selection process at each level if necessary.
pðxk jz1:k Þ Z 1 ¼ pðzk jxk Þ pðxk jxk1 Þpðxk1 jz1:k1 Þ dxk1 C
3.2. The explicit method
1 ¼ pðzk jxk Þpðxk jz1:k1 Þ C " N l 1 X pðzk jxk Þoil k1 dðxk fxl gikjk1 Þ ¼ C i¼1
ð76Þ
Explicit multires propagation: The unires particle propagation is given by
ð77Þ
xikþ1jk ¼ f ðxik Þ þ vik ;
XXX jXj 1 i¼1
m
pðxk jz1:k1 Þ ¼
~ ihk1 dðxk fxhjm gikjk1 Þ pðzk jxk Þo
¼
N X
oik1 dðxk xikjk1 Þ
i¼1 dTkjk1 ok1
#
¼
ð78Þ
ð86Þ
X
ðdjkjk1 ÞT ojk1
ð87Þ
j¼1
" N l 1 X pðzk jfxl gikjk1 Þoil k1 ¼ C i¼1
¼
L N=2 X
j¼1
ðTdðTT ðxjlhk xjlhkjk1 ÞÞÞT ojlhk1 , ð88Þ
N
þ
jXj 1 i¼1
m
~ ihk1 pðzk jfxhjm gikjk1 Þo
¼
where 2
#
sgnðfxhjm gikjk1 Þ Nl X
ð79Þ
oil k dðxk fxl gikjk1 Þ
i¼1 N hj
þ
XXX jXj 1 i¼1
m
~ ihk dðxk fxhjm gikjk1 Þ o
sgnðfxhjm gikjk1 Þ,
ð80Þ
where the updated weights are oil k ¼
1 pðzk jfxl gikjk1 Þoil k1 C
(81)
and ~ ihk o
(82)
The updated estimate vector can be derived from
x^ kjk ¼
oil k fxl gikjk1
þ
i¼1
sgnðfxhjm gikjk1 Þ.
N hj jX 1 þL X X j¼j 1 i¼1
6 6 xjlhkjk1 ¼ T6 6 4
f ðx1k1 Þ
3j
2
v1k1
3j
7 6 7 7 6 7 7 þ T6 ... 7 7 6 7 5 4 L 5 2L 2 f ðxk1 Þ vk1 .. .
ð89Þ
¼ T f ðTT xjlhk1 Þ þ vjlhk1
ð90Þ
¼ f T ðxjlhk1 Þ þ vjlhk1 .
ð91Þ
From Eqs. (86) to (87), a partition of N samples into N=2L blocks is used where L is the level number of the multires decomposition. Since it is impossible to analytically derive f T ðÞ in Eq. (91), we take another approach. ~ jlh , we can find the For the thresholded o k1 j corresponding x~ lhk1 . Therefore, we can just propagate the elements in x~ jk1 ¼ TT x~ jlhk1 which corre~ jk1 ¼ TT o ~ jlh . spond to the distinct elements in o k1
1 ~ ihk1 . ¼ pðzk jfxhjm gikjk1 Þo C
Nl X
ð85Þ
N=2L
sgnðfxhjm gikjk1 Þ
hj XX X
(84)
and the propagated density is described by
N hj
þ
vik Nð0; Qk Þ
m
~ ihk fxhjm gikjk1 o ð83Þ
~ jk1 , we can For the repeated elements of o propagate one representative element in x~ jk1 .The number of particle propagations is controlled by the thresholding operation, which can significantly reduce the computation in propagation while maintaining performance. The details of finding the corresponding elements in x~ jlhk1 given ~ jlh are described as follows. For an example o k1 of three levels (L ¼ 3), we can explicitly derive
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
the form of 2 3 h 6 3 6 h 6 6 2 6 h 6 6 6 0 6 T¼6 6 h 6 6 6 0 6 6 6 0 4
~ jlh ð2Þ ¼ 0()sum of 1st 4 and sum of 2nd 4 o k1 ~ jk1 are the same; elements of o j ~ lh ð3Þ ¼ 0()sum of 1st 2 and sum of 2nd 2 o k1 ~ jk1 are the same; elements of o j ~ lh ð4Þ ¼ 0()sum of 3rd and sum of 4th 2 o k1 ~ jk1 are the same; elements of o j ~ lh ð5Þ ¼ 0()1st and 2nd elements of o ~ jk1 o k1 are the same; ~ jlh ð6Þ ¼ 0()3rd and 4th elements of o ~ jk1 o k1 are the same; ~ jlh ð7Þ ¼ 0()5th and 6th elements of o ~ jk1 o k1 are the same; ~ jlh ð8Þ ¼ 0()7th and 8th elements of o ~ jk1 o k1 are the same.
T as h3
h3
h3
h3
h3
h3
h3
h3
h3
h3
h3
h3
h2
h2
h2
0
0 2
h2
0
0
0
h
h
0
0
0
0
0
0
h
h
0
0
0
0
0
0
h
h
0
h
7 h3 7 7 7 07 7 7 h2 7 7 7, 07 7 7 07 7 7 07 5
0 2
3
h3
0 0 0 0 0 0 h h pffiffiffi 2 h¼ . ð92Þ 2 Here is a list of correspondences between zero ~ jlh ~ jk1 : elements in o and the same elements in o k1
1397
A combination of the above correspondences could also result in a combination of the same ~ jlh ð2Þ ¼ 0 and element numbers. For instance, fo k1 j ~ lh ð7Þ ¼ 0g(){sum of 1st 4 and sum of 2nd 4 o k1 ~ jk1 are the same, and 5th and 6th elements of o
~ jlh ð1Þ ¼ 0()sum of 8 elements of o ~ jk1 is o k1 zeros;
original p (x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x unires propagated p (x1|0)
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
x unires updated p (x1|1)
p (x1|1)
0.03 0.02 0.01 0 -5
0
5
10
15
20
25
x Fig. 11. Traditional (unires) particle filtering with 1000 particles: (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1398
~ jk1 rm are the same} or fo ~ jlh ð2Þ ¼ elements of o k1
where wk Nð0; Rk Þ, we want to update the propagated density
~ jlh ð5Þ ¼ 0 and ojlh ð6Þ ¼ 0g(){sum of 1st 4 0; o k1 k1 ~ jk1 are the same, 1st and sum of 2nd 4 elements of o j ~ k1 are the same and 3rd and and 2nd elements of o ~ jk1 are the same}. 4th elements of o The propagation of the thresholded quantity is then given by x~ jkjk1 ¼ f ðTT x~ jlhk1 Þ þ vjk1
pðxk jz1:z1 Þ ¼
(93)
~ jk1 . dT ðxjk x~ jkjk1 Þo
l i ¼ Nðzk gðx~ ikjk1 Þ; Rk Þ;
~ i ¼ 1; . . . ; N,
(94)
Explicit multires update: Given the measurement model
pðxk jz1:z Þ ¼
zk ¼ gðxk Þ þ wk ,
where C ¼
(95)
N~ X 1 i i ~ k1 l dðxk x~ ikjk1 Þ, o C i¼1
PN~
~ ik1 i¼1 o
li .
original unires p (x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x
b
multires propagated p (x1|0), thresholded
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
15
20
x
c
multires updated p (x1|1), thresholded
p (x1|1)
0.04 0.02 0 -5
0
(97)
where N~ is the total number of distinct particles, which could be much smaller than N depending on the threshold. The density update can be carried out by
j¼1
a
(96)
with the explicit multires method. The likelihood functions are
N=2L
X
~ jk1 dT ðxjk x~ jkjk1 Þo
j¼1
which in turn defines the density propagation as pðxk jz1:k1 Þ ¼
L N=2 X
5
10 x
Fig. 12. Multires particle filtering with 562 particles (T s ¼ 0:0002): (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
(98)
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1399
Example 4. The nonlinear system and measurement models are given by
without the multiresolution approach, where
xk1 5xk1 þ xk ¼ þ vk ; 2 ð1 þ x2k1 Þ
x^ 0j0 ¼
vk Nð0; 1Þ
(99)
zk ¼ 10 tan1
k
10
þ wk ;
xi0 oi0
i¼1
and x
N X
wk Nð0; 1Þ
with xi0 px ðx0 Þ. Fig. 11 shows the propagation and update of px ðx0 Þ using the traditional particle filter approach (we call it the unires approach) with the number of particles being 1000. It can be seen that the particle cloud peak moves from 10 to 6 with a wider spread in (b) due to propagation and becomes more concentrated in (c) due to update. Figs. 12 and 14 show the propagation and update of px ðx0 Þ with the newly developed multires particle approach (only the explicit approach is used due to the paper length limit) where the threshold (T s ) is 0.0002 in Fig. 12, 0.001 in Fig. 13 and 0.002 in Fig. 14. The propagation-update cycle from t ¼ 0 to 1 can be applied to other time instants in any
(100)
with the initial nonlinear (bi-modal) state density being px ðx0 Þ ¼ 0:7 Nð10; 22 Þ þ 0:3 Nð20; 42 Þ.
(101)
In the following, we will show how px ðx0 Þ is propagated from t ¼ 0 to 1 using particles with and without the multiresolution approach. Then the propagated density will be updated by z1 ¼ 15 tan1 ðx^ 0j0 =10Þ þ w1 using particles with and
a
original unires p (x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x
b
multires propagated p (x1|0), thresholded
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
15
20
x
c
multires updated p(x1|1), thresholded
p (x1|1)
0.04 0.02 0 -5
0
5
10 x
Fig. 13. Multires particle filtering with 234 particles ðT s ¼ 0:001Þ: (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401
1400
a
original unires p(x0)
p (x0)
0.015 0.01 0.005 0 0
5
10
15
20
25
30
35
x
b
multires propagated p(x1|0), thresholded
p (x1|0)
0.015 0.01 0.005 0 0
5
10
15
20
25
x
c
multires updated p (x1|1), thresholded
p (x1|1)
0.04 0.02 0 -5
0
5
10
15
20
x Fig. 14. Multires particle filtering with 165 particles ðT s ¼ 0:002Þ: (a) px ðx0 Þ, (b) px ðx1j0 Þ, and (c) px ðx1j1 Þ.
Table 1 A comparison of computational complexity and estimation accuracy (500 runs) Approaches
Unires particle filter Multires particle filter ðT s ¼ 0:0002Þ Multires particle filter ðT s ¼ 0:001Þ Multires particle filter ðT s ¼ 0:002Þ
RMS estimation errors
Number of particles
0.0450 0.0451
1000 562
0.0443
234
0.0452
165
filtering application; therefore, a complete multires particle filtering algorithm is derived. A numerical comparison of estimation accuracy and computational complexity of both unires and multires approaches is given in Table 1. It can be clearly seen
that the proposed multires approach performs comparably (or even slightly better sometimes) to the unires approach, but is significantly more efficient! 4. Conclusions We have developed a framework for computationally efficient particle filtering using the spatial domain multiresolutional approach. Both implicit and explicit implementation methods for a multires particle filter are presented. Computational efficiency of a particle filter is achieved by drastically reducing the needed number of particles while maintaining a comparable estimation accuracy. References [1] N.J. Gordon, D.J. Slamond, A.F.M. Smith, Novel approach to nonlinear/non-Gaussian Bayesian state estimation, IEE Proc.-F 140 (2) (April 1993) 107–113.
ARTICLE IN PRESS L. Hong, D. Wicker / Signal Processing 87 (2007) 1384–1401 [2] R. Van der Merwe, N. de Freitas, A. Doucet, E.A. Wan, The unscented particle filter, Technical Report CUED/F-INFENG/TR380, Cambridge University Engineering Department, August 2000. [3] A. Doucet, N. de Freitas, N. Gordon, Sequential Monte Carlo Methods in Practice, Springer, New York, 2001. [4] N. Cui, L. Hong, J.R. Layne, A comparison of nonlinear filtering approaches with application to ground target tracking, Signal Processing 85 (8) (August 2005) 1469–1492. [5] A. Farina, B. Ristic, Tracking a ballistic target: comparison of several nonlinear filters, IEEE Trans. Aerospace Electronic Systems 38 (3) (July 2002) 854–867. [6] A. Doucet, N. Gordon, Efficient particle filters for tracking maneuvering targets in clutter, in: IEE Colloquium on Target Tracking: Algorithms and Applications (Ref. No. 1999/090, 1999/215), 1999, pp. 4/1–4/5.
1401
[7] D. Fox, KLD-sampling: adaptive particle filter, Adv. Neural Inf. Process. Syst. (NIPS) 14 (2001) 713–720. [8] L. Hong, Multiresolutional distributed filtering, IEEE Trans. Automat. Control 39 (4) (1994) 853–856. [9] L. Hong, Multiresolutional multiple-model target tracking, IEEE Trans. Aerospace Electronic Systems 30 (2) (April 1994) 518–524. [10] L. Hong, Multirate interacting multiple model filtering for target tracking using multirate models, IEEE Trans. Automat. Control 44 (7) (1999) 1326–1340. [11] L. Hong, N. Cui, M. Bakich, J.R. Layne, Multirate interacting multiple model particle filter with application to terrain-based ground target tracking, IEE Proc: Control Theory Appl. 153 (6) (2006) 721–731. [12] I. Daubechies, Orthonormal bases of compact supported wavelets, Commun. Pure Appl. Math. XLI (1988) 909–996.