Optimal and unbiased FIR filtering in discrete time state space with smoothing and predictive properties

We address p-shift finite impulse response optimal (OFIR) and unbiased (UFIR) algorithms for predictive filtering (p > 0), filtering (p = 0), and smoothing filtering (p < 0) at a discrete point n over N neighboring points. The algorithms were designed for linear time-invariant state-space signal models with white Gaussian noise. The OFIR filter self-determines the initial mean square state function by solving the discrete algebraic Riccati equation. The UFIR one represented both in the batch and iterative Kalman-like forms does not require the noise covariances and initial errors. An example of applications is given for smoothing and predictive filtering of a two-state polynomial model. Based upon this example, we show that exact optimality is redundant when N ≫ 1 and still a nice suboptimal estimate can fairly be provided with a UFIR filter at a much lower cost.


Introduction
There is a class of estimation problems requiring optimal filtering at a discrete-time current point n employing measurement on an averaging interval (horizon) of preceding or/and succeeding neighboring but not obligatorily nearest N points. To solve such problems, filtering is usually organized employing the finite impulse response (FIR) structures. Because the averaging interval can be placed with an arbitrary time shift p with respect to n, there can be recognized three kinds of p-shift FIR filters as shown in Figure 1, namely the p-step predictive filter (p > 0), filter (p = 0), and |p|-lag smoothing filter (p < 0).
Predictive FIR filtering is fundamental for discrete-time feedback systems and required in signal processing when measurement is temporary unavailable in the nearest past of p points. The one-step predictive filter known as the receding horizon filter has been put into the concept of the receding horizon (or model predictive) control [1,2]. For polynomial signals, an unbiased predictive FIR filter was proposed by Heinonen and Neuvo in [3]. Further, this *Correspondence: shmaliy@ugto.mx Electronics Department, DICIS, Universidad de Guanajuato, Salamanca, 36855, Gto., Mexico filter was investigated by many authors [4] and developed in state-space to p-step predictive filtering [5].
Smoothing FIR filtering is a key solution whenever denoising of signals is required with highest efficiency. Savitzky-Golay smoothing filter [6] is one of the most popular here. In recent decades, we meet a few new substantial results. Linear FIR smoothers were developed and used by Zhou and Wang in the FIR-median hybrid filters [7]. In state space, order-recursive FIR smoothers were proposed by Yuan and Stuller in [8]. Most recently, the general receding horizon FIR smoother theory has been developed in [9,10] and, for polynomial signals, the |p|-lag smoothing FIR filter theory addressed in [11].
It follows from the above-given short survey that the authors prefer solving the problems of filtering, smoothing, and prediction employing different algorithms. In [12,13], a universal scheme has been proposed for the pshift FIR estimators (filters, predictors, and smoothers). Still no universal solution has been addressed in state space for FIR filtering with smoothing and prediction properties.
In this article, we follow the approach developed in [12] and address universal p-shift optimal FIR (OFIR) http://asp.eurasipjournals.com/content/2012/1/163 … p-step predictive FIR filtering and unbiased FIR (UFIR) filters for predictive filtering (p > 0), filtering (p = 0), and smoothing filtering (p < 0) at a current point n of linear discrete-time-invariant statespace models with white noise. The rest of the article is organized as follows. In Section 'Signal model and problem formulation' , we describe the model and formulate the problem. The p-shift OFIR filter is derived in Section 'p-Shift OFIR filter with predictive and smoothing properties' . Here, we also find its gain and estimate the initial mean square state function. The UFIR filter is considered in detain in Section 'p-shift UFIR filter with predictive and smoothing properties' along with the estimation error. An application to the two-state model is given in Section 'Applications' and concluding remarks are drawn in Section 'Conclusion' .

Signal model and problem formulation
Consider a class of discrete time-invariant linear signal models represented in state space with the state and observation equations, respectively, where x n ∈ R K and y n ∈ R M are the state and observation vectors, respectively, Here, A ∈ R K×K , B ∈ R K×P , C ∈ R M×K , and D ∈ R M×M . The system noise vector w n ∈ R P and the measurement noise vector v n ∈ R M , respectively, are white Gaussian with zero mean components, E{w n } = 0 and E{v n } = 0. It is implied that w n and v n are mutually uncorrelated, E{w i v T j } = 0 for all i and j, and have the covariances, respectively, The problem now formulates as follows. Given the model (1) and (2), we would like to derive a p-shift OFIR filter covering the problems of predictive filtering (p > 0), filtering (p = 0), and smoothing filtering (p < 0) as shown in Figure 1. We also wish to find its unbiased version, represent it in the iterative Kalman-like form, and investigate errors based on a typical example.

p-Shift OFIR filter with predictive and smoothing properties
A p-shift OFIR filter can be derived following Figure 1, if to represent (1) and (2) on a horizon of N points, similarly http://asp.eurasipjournals.com/content/2012/1/163 to [11], with recursively computed forward-in-time solutions [14] as follows, respectively, where X n,m ∈ R KN , Y n,m ∈ R MN , W n,m ∈ R PN , and V n,m ∈ R MN are given by, respectively, The matrices where we have assignedC In this model, the initial state x m−p is supposed to be known exactly and w m−p is thus allowed to have zero components.

Optimal gain
One can now assign the gain matrix H(p) H(p, n, m) ∈ R K×MN implementing the convolution principle and find the filtering estimate a of x n as For H(p) to be optimal in the minimize mean square error (MSE) sense, the following cost function needs to be minimized, where E(x) means an average of x. The minimization can be provided employing the orthogonality condition [14] in the form of [12], to produce the optimal gain matrixĤ(p). In doing so, one needs substituting x n with the first vector row b in (9); that is, whereB N−1 is the first vector row in (16). Substituting (24) to (23) and supposing that the initial state and measurement noise are mutually uncorrelated for all p, we provide the averaging in (23) and arrive at the optimal gain matrix the mean square initial state is specified by By multiplying R m−p in (25) from the left-hand side with the identity matrix (32) where, by n − m = N − 1, the unbiased gain attains two equivalent forms of Note thatH(p) satisfies the unbiasedness condition and has an important applied property: it does not depend on noise and initial errors, although it is p-and Ndependent. As shown in Appendix, matrix Z m−p representing in (32) the mean square initial state R m−p on an averaging interval of N points can optimally be estimated by solving the discrete algebraic Riccati equation (DARE) whose analytic solution can be found following [15]. We notice that this equation also serves for filtering out all of the noise components [12].

Optimal filtering estimate
Determined Z m−p , by (36), the p-shift OFIR filtering estimatex n|n−p can now be generalized as follows. Given (1) and (2) with uncorrelated zero-mean white noise vectors w n and v n . Then p-step OFIR predictive filtering (p > 0), filtering (p = 0), and |p|-lag smoothing filtering (p < 0) can be provided at n employing data taken from m − p to n − p byx where Y n−p,m−p is the data vector (12), C i is given by (17), andZ w ,Z w , andZ v are specified for (25). The algorithm should be applied to any N K, in order to avoid problems with singularities. Note that K is typically not larger in state space modeling.

p-shift UFIR filter with predictive and smoothing properties
There are at least two cases when exact optimality is redundant and OFIR can fairly be substituted with UFIR at much lower price to produce still a nice near optimal estimate [12]. In fact, if Z m−n substantially dominatesZ w , Z w , andZ v in the order of magnitudes for all p, we havê H(p) ∼ =H(p). The same effect is achieved with N 1. Thus, the UFIR filter should also be generalized. Given (1) and (2) with uncorrelated zero-mean white noise components w n and v n . Then p-step UFIR predictive filtering (p > 0), filtering (p = 0), and |p|-lag smoothing filtering (p < 0) can be provided at n employing data taken from m − p to n − p bȳ Note that both (40) and (41) follow from (38) and (39) straightforwardly if to refer to the unbiasedness condition (35) and neglectZ w ,Z w , andZ v .

Kalman-like UFIR filtering algorithm
Noticing that the UFIR filter (41) ignoring the noise covariances and initial error is highly attractive for engineering applications, one also notes that the computational problem may arise in its batch form when N 1. To circumvent this problem, a fast iterative Kalman-like form has been addressed in [13] for filters, predictors, and smoothers. If to introduce a time shift p tox n+p|n stated by Theorem 2 in [13] for time-invariant models and take into consideration that initial F l is time-invariant, then the iterative Kalman-like form of (41) appears as follows: http://asp.eurasipjournals.com/content/2012/1/163 where K l K l (p) = A p F l C T is the bias correction gain, m = n − N + 1, s = α − 1, and an iterative variable l ranges from α = max(m + K, m + 2, m + 2 − p) to n. The true estimate corresponds to l = n.
As well as in the case of UFIR estimator [13], the gain K l in (46) also does not depend on noise and initial errors. In this algorithm, we have two batch forms, (43) and (44), which can be computed fast for typically small K. To avoid singularities, the computation starts with l = α and finishes at l = n. This last value is used as true and the procedure repeated iteratively for each new measurement. The iterative p-shift Kalman-like algorithm (42)-(46) is listed in Table 1 in a convenient computation form.

Estimation error
Although the estimation error is not involved to the algorithm (42)-(46) that is its extremely remarkable property, the MSE may be required to characterize the filter performance.
The MSE in the p-shift FIR filtering estimate can be evaluated by the matrix Substitutingx l|l−p with (46), assigningx l x l|l−p and ε l = x l −x l , and employing (1) and (2), we first write Usex l|l−p as the output when l = n As a next step, it needs to express x l−p via x l−1 . That can be done if we write (1) forward and backward in time for different p and provide the transformations in order to have finally where By substituting (49) with (50) to (48), taking into consideration that E{ε l−1 β T l } and E{β l ε T l−1 } have zero components, and providing the averaging, we finally come up with whereR = E{w l β T l },R = E{β l w T l }, andR = E{β l β T l } are given with, respectively, By (51), the estimation error can thus be computed iteratively along with the estimate (46). As can be seen, P n inherently diminishes in smoothing filtering (p < 0), bȳ R andR. It rises with higher rate in predictive filtering http://asp.eurasipjournals.com/content/2012/1/163 (p > 0) owing to the effect ofP. In the case of filtering (p = 0), P n is computed by In all of the cases, P n becomes zero if the model is deterministic and the filter order is exactly that of a system.
Note that P n computed in such a way ranges upper the true value due to the accumulating effect caused by the noise covariances. Alternatively, P n can be well bounded with the error bound (EB) specified in [16] in the threesigma sense via the noise power gain (NPG) as to characterize the noise standard deviation in the vto−g filter channel via measurement of the kth state in the presence of white noise having the variance σ 2 k . The NPG coefficient g k (vg) g k(vg) (N, p) is defined here as a components of the NPG matrix K k K k (N, p) (11) where the thinned gain H k H (p) k ∈ R K×N is composed with the Kth columns ofH(p) given by (42) starting with the kth one as To avoid the computational problem with N 1, the NPG K k can be computed iteratively [12] as by changing an iterative variable j from j = γ K, to N − 1. The initial value K k(γ −1) is provided by (58), if to substitute N with γ , and the true K k is taken when l = N − 1.

Applications
A comparison of errors in the FIR and Kalman filters has been provided in many articles [2,9,10,12,13]. Much lesser attention has been paid to the trade-off between the OFIR and UFIR filter outputs. To investigate errors in the proposed p-shift OFIR and UFIR filters and thereby learn their facilities, below we exploit a two-state model represented with (1) and (2) (7) and (8) of zero mean noise components, w n and v n , are allowed to be , respectively. Measurement has been provided in the presence of the zero-mean noise v n uniformly distributed from −2 to 2 with the variance σ v = 2/ √ 3. Both the OFIR algorithm and UFIR one ( Table 1) have been tested and the filtering errors evaluated in the first state at a current point n for different p and fixed N. Errors were bounded with EB β(p) β 1(11) (N, p) calculated by (56).

Errors in predictive FIR filtering
Supposing that the estimate is required at n = 50 and assuming that measurement may not be available in the nearest past points (as it sometimes occurs in wireless systems), we let 0 p 30 and find the predictive filtering estimate for N = 10 and N = 20. Figure 2 sketches errors in OFIR and UFIR estimates as functions of p.
Here, we also show the bounds ±β(p) for each N as functions of p. Note that for the model in question, EB can also be calculated via NPG g 1 (N, p) found in [5] as β 1(11) (N, p) Observing Figure 2, one infers that the estimates are closely related and that the errors range well within EBs stretched by growing p. Inherently, the prediction error is reduced by increasing N that can be seen by comparing Figure 2a,b.

Errors in smoothing FIR filtering
In the second experiment, we change p from zero to −N + 1 and evaluate errors in the smoothing filters at n = 30. Figure 3 illustrates the results for N = 10 ( Figure 3a) and N = 20 (Figure 3b).
The first conclusion that can be made is that errors reach a minimum at a center of the averaging horizon with p = −N/2. This should not be surprising, since the ramp impulse response associated with linear models reduces at this point to the uniform one associated with simple averaging [11] producing noise minimum among all other filters [6]. It can also be seen that the difference between the optimal and unbiased estimates exists but it is not http://asp.eurasipjournals.com/content/2012/1/163

Conclusion
In this article, the p-shift OFIR and UFIR algorithms with the properties of predictive filtering (p > 0), filtering (p = 0), and smoothing filtering (p < 0) have been addressed for linear discrete time-invariant state-space models. The OFIR filter is shown to self-determine the mean square initial state function by solving the DARE. The UFIR filter represented both in the batch and iterative Kalman-like forms ignores covariances and initial errors, unlike the Kalman filter. As an example of applications, we have exploited the two-state polynomial model and investigated errors in the OFIR and UFIR filters. Based upon, we have confirmed the statement made earlier for FIR filters, predictors and smoothers: the UFIR estimate converges to the OFIR one by increasing N and the estimation errors are well bounded with EB. That means that exact optimality may be redundant with N 1 and still a nice suboptimal estimate can be provided with UFIR filters at much lower cost. An importance of the OFIR and UFIR filtering algorithms proposed resides in the fact that they are both general for linear discrete time-invariant state-space models. The algorithms virtually generalize the well-known Savitzky-Golay solution for smoothing and predictive filtering in state-space. However, unlike the latter, both OFIR and UFIR filters have the convolution-based forms more familiar for electronics engineers. Moreover, the convolution computation can efficiently be provided in the frequency domain that is its another benefit. Finally, engineers should certainly appreciate the iterative Kalman-like algorithm. Paying attention to these advantages, our current investigations are focused on several applied problems associated with signal and image processing.
Endnotes ax n|k is the filtering estimate at n via measurement from the past to k;x n|k is optimal andx n|k unbiased. b The case of filtering out all of the noise components is considered in [12].