AKALMANLIKEFIRESTIMATORIGNORINGNOISEANDINITIALCONDITIONS
Yuriy S. Shmaliy
Electronics Department, Guanajuato UniversityCtra. SalamancaValle, 3.5+1.8km, PaloBlanco, 36855, Salamanca, Mexicophone: + 52 (464) 6470195, email: shmaliy@ugto.mx,web: www.ingenierias.ugto.mx
ABSTRACT
A
p
shift ﬁnite impulse response (FIR) unbiased estimator(UE) is addressed for linear discrete timevarying ﬁltering(
p
=
0),
p
step prediction (
p
>
0), and
p
lag smoothing(
p
<
0) of signal models in state space with no requirements for initial conditions and zero mean noise. A solutionis found in a batch form and represented in a computationally efﬁcient iterative Kalmanlike one. It is shown that theKalmanlike FIR UE is able to overperform the Kalman ﬁlterif the noise covariances and initial conditions are not knownexactly, noise is not white, and both the system and measurement noise components need to be ﬁltered out. Otherwise,the errors are similar.
1. INTRODUCTION
For such unsuited applications of the Kalman ﬁlter (KF) [1]as estimation of nonlinear models, under unknown initialconditions, and in the presence of nonwhite or multiplicative noise sources, the Kalmanlike one is often designed tosave the recursive structure, while connecting the algorithmcomponents with the model in different ways. Because therecan be found an inﬁnity of Kalmanlike solutions depending on applications, we meet a number of propositions suggesting some new qualities while saving (or not deterioratingsubstantially) the advantages of KF: fast computation and accuracy.Cox in [2] and others have derived the extended KF(EKF) for nonlinear models by a linearization of the statespace equations. Referring to the fact that EKF can give particularly poor performance when the model is highly nonlinear [3], Julier and Uhlmann employed in [4] the unscentedtransform and proposed the unscented KF (UKF). Both EKFand UKF have then been used extensively and the formerwas developed in [5] to the invariant EKF for nonlinearsystems possessing symmetries (or invariances). For highdimensional systems, the ensemble KF was proposed byEvensen in [6] and, for systems with sparse matrices, the fastKF applied by Lange in [7]. Applications has also found therobust Kalmantype ﬁlter designed by Masreliez [8] [9] forlinear statespace relations with nonGaussian noise referredto as heavy tailed noise or Gaussian one mixed with outliers.UsefulKalmanlikealgorithmscanalsobefoundinworksbyNahi [10], Basseville
et al
. [11], Baccarelli and Cusani [12],AitElFquih and Desbouvries [13], Carmi
et al
. [14], Stefanatos and Katsaggelos [15], and the list can be extended.In spite of great efforts in extending the applications andimproving the performance of KF, its structure still remainsrecursive thus with the inﬁnite impulse response (IIR). Investigating in [3] both the IIR and ﬁnite impulse response(FIR) ﬁlters, Jazwinski resumed that the limited memory ﬁlter (having FIR) appears to be more robust against the unbounded perturbation in the system. Referring to [3], optimal FIR ﬁltering has been developed by W. H. Kwon
et al
.in [16]. There were also proposed several Kalmanlike FIRestimators by Kwon
et al
. in [17], Han
et al
. in [18], andShmaliy in [19]. A distinctive feature of such algorithms isthat white Gaussian noise in the convolutionbased estimateobtained over
N
past measured points is reduced as a reciprocal of
N
[20] disregarding the model [19]. Moreover, the unbiased and optimal FIR estimates typically become stronglyconsistent if
N
occurs to be large [21] or the mean square initial state function dominates the noise covariance functionsin the order of magnitudes [19]. It is also known that the optimal horizon
N
opt
makes the FIR estimate (optimal or unbiased) similar or even better than the Kalman one [16,19–23].Owing to the exciting engineering features of theKalmanlike FIR algorithms uniting advantages of KF andinherent properties of FIR structures such as the boundedinput/bounded output (BIBO) stability as well as better robustness against temporary model uncertainties, non Gaussian noise, and roundoff errors, it may be expected that theFIR unbiased estimator (UE) ignoring noise and initial conditionswillserveefﬁciently insteadofoptimal ﬁltersinmanyapplications.
2. SIGNAL MODEL
Consider a class of discrete timevarying (TV) linear statespace models represented with the state and observationequations, respectively,
x
n
=
A
n
x
n
−
1
+
B
n
w
n
,
(1)
y
n
=
C
n
x
n
+
D
n
v
n
,
(2)where
x
n
∈
ℜ
K
and
y
n
∈
ℜ
M
are the state and observationvectors, respectively. Here,
A
n
∈
ℜ
K
×
K
,
B
n
∈
ℜ
K
×
P
,
C
n
∈
ℜ
M
×
K
, and
D
n
∈
ℜ
M
×
M
. The vectors
w
n
∈
ℜ
P
and
v
n
∈
ℜ
M
are zero mean,
E
{
w
n
}
=
0
and
E
{
v
n
}
=
0
. It is impliedthat
w
n
and
v
n
are mutually uncorrelated and independentprocesses,
E
{
w
i
v
T j
}
=
0
, having arbitrary distributions andknown covariances
Q
w
(
i
,
j
) =
E
{
w
i
w
T j
}
,
(3)
Q
v
(
i
,
j
) =
E
{
v
i
v
T j
}
,
(4)for all
i
and
j
, to mean that
w
n
and
v
n
should not obligatorilybe Gaussian and deltacorrelated.Following the strategies of the recursive KF [1] and iterative Kalmanlike FIR unbiased ﬁlter (UF) [19], the TV
19th European Signal Processing Conference (EUSIPCO 2011)Barcelona, Spain, August 29  September 2, 2011
© EURASIP, 2011  ISSN 20761465985
A
n
Gain
Actual signal
Time
A
n
–
N
+2
A
n
–
N
+2
Gain
Averaging horizon of
N
points
A
n
E s t i m a t e
Recursive s
trategy of the
KF
Iterative s
trategy of the Kalman

like FIR UF
n N+

1
n
opt
Figure 1: Strategies of the recursive KF and iterativeKalmanlike FIR UF algorithms.estimates of
x
n
can be obtained as shown in Fig. 1. Here,KF starts at some initial point
n
−
N
+
1, where
N
2, withknown initial conditions and recursively produces estimatesat each subsequent point up to
n
. The estimate is formed intwo steps. First, the system matrix
A
n
−
N
+
2
makes a projection from
n
−
N
+
1 to
n
−
N
+
2 and then the Kalman gainadjusts the result to be the estimate at
n
−
N
+
2. The procedure repeats recursively, provided the covariances of whitenoise sequences.The Kalmanlike FIR UF does not require the noise performance and initial conditions [24], but needs an optimalaveraging horizon
N
opt
[25] in order for the estimate to havethe minimum mean square error (MSE) and be consistent tothe Kalman one [19,22]. This ﬁlter starts with any unknownvalue at
n
−
N
+
1 and iteratively produces estimates at subsequent points in two steps similarly to KF, although the truevalue is taken only at
n
when
N
=
N
opt
. The algorithm operates in any noise environment that makes it highly attractivefor engineering applications.
3. TIMEVARYING BATCH UNBIASED FIRESTIMATOR
In order to ﬁnd the FIR UE, (1) and (2) can be extended on ahorizon of
N
points from
m
=
n
−
N
+
1 to
n
following [26]and similarly to [21] as, respectively,
X
n
,
m
=
A
n
,
m
x
m
+
B
n
,
m
W
n
,
m
,
(5)
Y
n
,
m
=
C
n
,
m
x
m
+
G
n
,
m
W
n
,
m
+
D
n
,
m
V
n
,
m
,
(6)where
X
n
,
m
∈
ℜ
KN
,
Y
n
,
m
∈
ℜ
MN
,
W
n
,
m
∈
ℜ
PN
, and
V
n
,
m
∈
ℜ
MN
are speciﬁed by, respectively,
X
n
,
m
=
x
T n
x
T n
−
1
...
x
T m
T
,
(7)
Y
n
,
m
=
y
T n
y
T n
−
1
...
y
T m
T
,
(8)
W
n
,
m
=
w
T n
w
T n
−
1
...
w
T m
T
,
(9)
V
n
,
m
=
v
T n
v
T n
−
1
...
v
T m
T
,
(10)and
A
n
,
m
∈
ℜ
KN
×
K
,
C
n
,
m
∈
ℜ
MN
×
K
,
G
n
,
m
∈
ℜ
MN
×
PN
, and
D
n
,
m
∈
ℜ
MN
×
MN
are given with, respectively,
A
n
,
m
=
A
m
+
1
T
n
A
m
+
1
T
n
−
1
...
A
T m
+
1
I
T
,
(11)
C
n
,
m
=
C
n
A
m
+
1
n
C
n
−
1
A
m
+
1
n
−
1
...
C
m
+
1
A
m
+
1
C
m
,
(12)
G
n
,
m
=
¯
C
n
,
m
B
n
,
m
,
(13)
D
n
,
m
=
diag
D
n
D
n
−
1
...
D
m
N
,
(14)where we assigned
A
n
−
gn
−
h
=
g
∏
i
=
h
A
n
−
i
,¯
C
n
,
m
=
diag
C
n
C
n
−
1
...
C
m
N
, and performed
B
n
,
m
∈
ℜ
KN
×
PN
as
B
n
,
m
=
B
n
A
n
B
n
−
1
...
A
m
+
2
n
B
m
+
1
A
m
+
1
n
B
m
0 B
n
−
1
...
A
m
+
2
n
−
1
B
m
+
2
A
m
+
1
n
−
1
B
m
+
1
...............
0 0
...
B
m
+
1
A
m
+
1
B
m
0 0
...
0 B
m
,
(15)The model, (5) and (6), suggests that the state equation atthe initial point
m
is
x
m
=
x
m
+
B
n
w
m
that for
B
n
specialized with (15) can uniquely be satisﬁed with
w
m
zerovalued.The initial state
x
m
must thus be known
a priori
or estimatedoptimally
a posteriori
as shown in [19].By the convolution, the estimate
1
˜
x
n
+
p

n
of
x
n
can nowbe obtained if we assign a
K
×
MN
gain matrix
H
n
,
m
(
p
)
andclaim that˜
x
n
+
p

n
=
H
n
,
m
(
p
)
Y
n
,
m
(16a)
=
H
n
,
m
(
p
)(
C
n
,
m
x
m
+
G
n
,
m
W
n
,
m
+
D
n
,
m
V
n
,
m
)
.
(16b)The estimate (16a) will be unbiased if and only if thefollowing unbiasedness condition is satisﬁed
E
{
˜
x
n
+
p

n
}
=
E
{
x
n
+
p
}
,
(17)where
E
means averaging of the succeeding relation.Averaging in (16b), by (17), means removing thezero mean noise matrices that gives us
E
{
˜
x
n
+
p

n
}
=
¯
H
n
,
m
(
p
)
C
n
,
m
x
m
, where¯
H
n
,
m
(
p
)
is the FIR UE gain. In turn,
E
{
x
n
}
can be substituted with the ﬁrst vector row of (5) byremoving noise as
E
{
x
n
}
=
A
m
+
1
n
x
m
. Since
n
can be arbitrary, one can substitute it with
n
+
p
and write
E
{
x
n
+
p
}
=
A
m
+
1
n
+
p
x
m
.
(18)Equating
E
{
˜
x
n
+
p

n
}
to (18) leads to the unbiasednessconstraint for TV models
A
m
+
1
n
+
p
=
¯
H
n
,
m
(
p
)
C
n
,
m
.
(19)
1
˜
x
n
+
p

n
is an estimate at
n
+
p
via measurement from the past to
n
; ˆ
x
n
+
p

n
mean optimal and ¯
x
n
+
p

n
unbiased.
986
If we further multiply (19) from the right hand sides withthe identity matrix
(
C
T n
,
m
C
n
,
m
)
−
1
C
T n
,
m
C
n
,
m
and then remove
C
n
,
m
from both sides, we go to¯
H
n
,
m
(
p
) =
A
m
+
1
n
+
p
(
C
T n
,
m
C
n
,
m
)
−
1
C
T n
,
m
(20)representing the gain of the TV FIR UE. It can easily be veriﬁed that (20) becomes that derived in [19] for timeinvariant(TI) models.Provided (20), the TV batch FIR UE is speciﬁed by thefollowing theorem, which proof belongs to (5)(20).
Theorem 1
Given (1) and (2) with zero mean mutually uncorrelated and independent
w
n
and
v
n
having arbitrary distributions and known covariance functions. Then, ﬁltering
(
p
=
0
)
, plag smoothing (p
<
0
), and pstep prediction (p
>
0
) are provided at n
+
p using data taken fromm
=
n
−
N
+
1
to n by the batch FIR UE as
¯
x
n
+
p

n
=
¯
H
n
,
m
(
p
)
Y
n
,
m
(21a)
=
A
m
+
1
n
+
p
(
C
T n
,
m
C
n
,
m
)
−
1
C
T n
,
m
Y
n
,
m
,
(21b)
where
C
n
,
m
is given by (12) and
Y
n
,
m
is the data vector (8).
4. TIMEVARYING KALMANLIKE ESTIMATOR
Although theorem 1 establishes an exact convolutionbasedrule to estimate unbiasedly the TV state at
n
as shown in Fig.1, the computational problem arises when
N
1 owing tolarge dimensions of all of the matrices and vectors. For fastcomputation, the batch FIR UE (21b) can be represented inan iterative Kalmanlike form stated by the following theorem, which proof is similar to that given in [19].
Theorem 2
Given the batch FIR UE (theorem 1). Then itsiterative Kalmanlike form is the following:
¯
x
l
+
p

l
=
A
l
+
p
¯
x
l
+
p
−
1

l
−
1
+
A
l
+
p
Υ
−
1
l
(
p
)
F
l
C
T l
×
[
y
l
−
C
l
Υ
l
(
p
)
¯
x
l
+
p
−
1

l
−
1
]
,
(22)
in which
Υ
l
(
p
) =
A
l
−
p

l
,
p
−
1
(
smoothing
)
A
l
,
p
=
0
(
ﬁltering
)
I
,
p
=
1
(
prediction
)
p
−
1
∏
i
=
1
A
−
1
l
−
i
,
p
>
1
(
prediction
)
,
F
l
= [
C
T l
C
l
+(
A
l
F
l
−
1
A
T l
)
−
1
]
−
1
,
(23)¯
x
s
+
p

s
=
A
m
+
1
s
+
p
PC
T s
,
m
Y
s
,
m
,
(24)
F
s
=
A
m
+
1
s
P
A
m
+
1
T
s
,
(25)
P
= (
C
T s
,
m
C
s
,
m
)
−
1
,
(26)
where s
=
m
+
K
−
1
, m
=
n
−
N
+
1
, and an iterative variablel ranges from m
+
K to n, because
C
T l
,
m
C
l
,
m
is singular withl
<
m
+
K. The true estimate corresponds to l
=
n.
As can be seen, (22) is the Kalman estimate, in which
A
l
+
p
Υ
−
1
l
(
p
)
F
l
C
T l
plays the role of the Kalman gain that,however, does not depend on noise and initial conditions.The algorithm has two batch forms, (24) and (25), which canbe computed fast for small
K
.Table 1: FullHorizon TV KalmanLike FIR UE AlgorithmStageGiven:
K
,
p
,
n
K
Set:
Υ
n
(
p
)
by (23)
P
= (
C
T K
−
1
,
0
C
K
−
1
,
0
)
−
1
F
K
−
1
=
A
1
K
−
1
P
A
1
T
K
−
1
¯
x
K
+
p
−
1

K
−
1
=
A
1
K
+
p
−
1
PC
T K
−
1
,
0
Y
K
−
1
,
0
Update:
F
n
= [
C
T n
C
n
+(
A
n
F
n
−
1
A
T n
)
−
1
]
−
1
¯
x
n
+
p

n
=
A
n
+
p
¯
x
n
+
p
−
1

n
−
1
+
A
n
+
p
Υ
−
1
n
(
p
)
F
n
C
T n
×
[
y
n
−
C
n
Υ
n
(
p
)
¯
x
n
+
p
−
1

n
−
1
]
4.1 FullHorizonTimeVaryingKalmanLikeEstimator
In special cases when noise is nonstationary or both the system and measurement noise components need to be ﬁlteredout, all the data available should be processed. By letting
N
=
n
+
1 and
l
=
n
K
in (22)–(26), the relevant fullhorizon algorithm becomes as shown in Table 1 and, for TImodels, simpliﬁes to that proposed in [19]. As can be seen,the algorithm (Table 1) requires only
K
and
p
, thus has extremely strong engineering features.
4.2 Error Bound
Provided¯
H
n
,
m
(
p
)
, the estimate error bound can be ascertained via the noise power gain (NPG) in the threesigmasense as follows:
EB
k
(
vg
)
(
n
,
N
,
p
) =
3
σ
k
K
1
/
2
k
(
vg
)
(
n
,
N
,
p
)
(27)where
σ
k
is the noise variance of the measurement of the
k
thstateand
K
k
(
vg
)
(
n
,
N
,
p
)
isthe
(
vg
)
thcomponentofthesquareNPG matrix
K
k
K
k
(
n
,
N
,
p
)
specialized as
K
k
=
H
k
H
T k
.Here the thinned
K
×
N
gain
H
k
¯
H
n
,
m
(
p
)
k
is composedby
K
th columns of ¯
H
n
,
m
(
p
)
starting with the
k
th one.
5. EXAMPLES OF APPLICATIONS
Below, we provide ﬁltering with
p
=
0 and prediction with
p
>
0 of the twostate polynomial model, (1) and (2), speciﬁed with
B
n
=
I
,
D
n
=
I
,
C
n
= [
1 0
]
, and
A
n
=
1
(
1
+
d
n
)
τ
0 1
,
(28)where
d
n
temporary takes different values. Such a situation occurs in oscillators undergoing temporary frequency“jumps” or in moving vehiculars with velocity “jumps”. ForTI ﬁltering,
d
n
represents uncertainty and, in the TV case,
d
n
is supposed to be known exactly. We mostly allow noise tobe white Gaussian, noticing that the relevant investigationsfor the uniformly distributed and highly intensive sawtoothnoise were provided in [20–22,24,29,30].
987
REFERENCES
[1] R. E. Kalman, “A new approach to linear ﬁltering andprediction problems,”
J. Basic Engineer.
, vol. 82, pp.35–45, Mar. 1960.[2] H. Cox, “On the estimation of state variables and parameters fornoisy dynamic systems,”
IEEE Trans. Autom. Contr.
, vol. 9, pp. 5–12, Jan. 1964.[3] A. H. Jazwinski,
Stochastic Processes and FilteringTheory
, New York: Academic Press, 1970.[4] S. J. Julier and J. K. Uhlmann, “A new extension of theKalmanﬁltertononlinearsystems”,
Proc.ofSPIE
,vol. 3068, pp. 182193, 1997.[5] S. Bonnabel, Ph. Martin and E. Sala¨un, “InvariantExtended Kalman Filter: theory and application to avelocityaided attitude estimation problem”, in
Proc.48th IEEE Conf. Decision Contr.
, pp. 1297–1304,2009.[6] G. Evensen, “Sequential data assimilation with nonlinear quasigeostrophic model using Monte Carlomethods to forecast error statistics,”
Journal of Geo physical Research
, vol. 99, pp. 143–162, May 1994.[7] A. A. Lange, “Simultaneous statistical calibration of the GPS signal delay measurements with related meteorological data”,
Physics and Chemistry of the Earth,Part A: Solid Earth and Geodesy
, vol. 26, pp. 471–473, 2001.[8] C. J. Masreliez, “Approximate nonGaussian ﬁltering with linear state and observation relations,”
IEEE Trans. Autom. Contr.
, vol. AC20, pp. 107–110, Feb1975.[9] R. D. Martin and C.J. Masreliez, “Robust estimationvia stochastic approximation,”
IEEE Trans. Inform.Theory
, vol. IT21, pp. 263–271, May 1975.[10] N. E. Nahi, “Optimal recursive estimation with uncertain observation,”
IEEE Trans. Inf. Theory
, vol. IT15,pp. 457–462. Jul. 1969.[11] M. Basseville, A. Benveniste, K. C. Chou, S. A.Golden, R. Nikoukhah, and A. S. Willsky, “Modeling and estimation of multiresolution stochastic processes,”
IEEE Trans. Inf. Theory
, vol. 38, pp. 766784,Mar. 1992.[12] E. Baccarelli and R. Cusani, “Recursive Kalmantypeoptimal estimation and detection of hidden Markovchains,”
Signal Process.
, vol. 51, pp. 55–64. May1996.[13] B. AitElFquih and F. Desbouvries, “Kalman Filtering in Triplet Markov Chains,”
IEEE Trans. SignalProcess.
, vol. 54, pp. 2957–2963, Aug. 2006.[14] A. Carmi, P. Gurﬁl, and D. Kanevsky, “Methods forsparse signal recovery using Kalman ﬁltering withembedded pseudomeasurement norms and quasinorms”,
IEEE Trans. on Signal Process.
, vol. 58, pp.2405–2409, Apr. 2010.[15] S. Stefanatos and A. K. Katsaggelos, “Joint data detection and channel tracking for OFDM systems withphase noise,”
IEEE Trans. Signal Process.
, vol. 56, pp.4230–4243, Sep. 2008.[16] W. H. Kwon and S. Han,
Receding horizon control:model predictive control for state models
, London:Springer, 2005.[17] W. H. Kwon, P. S. Kim, and P. Park, “A receding horizon Kalman FIR ﬁlter for discrete timeinvariant systems,”
IEEE Trans. Autom. Contr.
vol. 44, pp. 1787–1791, Sep. 1999.[18] S. H. Han, W. H. Kwon, and P. S. Kim, “Quasideadbeat minimax ﬁlters for deterministic statespacemodels,”
IEEE Trans. Autom. Contr.
, vol. 47, pp.1904–1908, Nov. 2002.[19] Y. S. Shmaliy, “Linear optimal FIR estimation of discrete timeinvariant statespace models,”
IEEE Trans.Signal Process.
, vol. 58, pp. 3086–3096, Jun. 2010.[20] Y. S. Shmaliy, “An unbiased FIR ﬁlter for TIE modelof a local clock in applications to GPSbased timekeeping,”
IEEE Trans. on Ultrason., Ferroel. and Freq. Contr
., vol. 53, pp. 862–870, May 2006.[21] Y. S. Shmaliy, “Optimal gains of FIR estimators for aclass of discretetime statespace models,”
IEEE Signal Process. Letters
, vol. 15, pp. 517–520, 2008.[22] P. S. Kim, “An alternative FIR ﬁlter for state estimationindiscretetimesystems,”
DigitalSignalProcess.,
vol. 20, pp. 935–943, May 2010.[23] W. H. Kwon, P. S. Kim, and S. H. Han, “A receding horizon unbiased FIR ﬁlter for discretetime statespace models,”
Automatica
, vol. 38, pp. 545–551,Mar. 2002.[24] Y.S.Shmaliy, “UnbiasedFIRﬁlteringofdiscretetimepolynomial statespace models,”
IEEE Trans. SignalProcess.
, vol. 57, pp. 1241–1249, Apr. 2009.[25] Y. S. Shmaliy, J. Mu˜nozDiaz, and L. ArceoMiquel,“Optimal horizons for a oneparameter family of unbiased FIR ﬁlters,”
Digital Signal Process.,
vol. 18, pp.739–750, Sep. 2008.[26] H. Stark and J. W. Woods,
Probability, Random Processes, and Estimation Theory for Engineers
, 2
nd
ed.,Upper Saddle River, NJ: Prentice Hall, 1994.[27] Y. S. Shmaliy, O. IbarraManzano, L. ArceoMiquel,and J. Mu˜nozDiaz, “A thinning algorithm for GPSbased unbiased FIR estimation of a clock TIE model,”
Measurement,
vol. 41, pp. 538–550, Jun. 2008.[28] M. Blum, “On the mean square noise power of an optimum linear discrete ﬁlter operating on polynomialplus white noise input,”
IRE Trans. Information Theory
, vol. 3, pp. 225–231, Dec. 1957.[29] Y. S. Shmaliy and O. IbarraManzano, “Optimal FIRﬁltering of the clock time errors,”
Metrologia
, vol. 45,pp. 571–576, Sep. 2008.[30] Y. S. Shmaliy, “A simple optimally unbiased MA ﬁlter for timekeeping,”
IEEE Trans. Ultrason. Ferroel.Freq. Control
, vol. 49, pp. 789–797, Jun. 2002.[31] Y. S. Shmaliy, “Linear unbiased prediction of clock errors,”
IEEE Trans. Ultrason. Ferroel. Freq. Control
,vol. 56, pp. 2027–2029, Sep. 2009.[32] G. H. Golub and G. F. van Loan,
Matrix Computations,
3
rd
Ed., Baltimore: The John Hopkins Univ.Press, 1996.
989