State-Space Representation of Forward and Backward Filters

Question

In [1], the author shows an efficient way of implementing the forward and backward filter using matrices. One can also implement this using filtfilt command in MATLAB. However, I am interested in implementing this using matrices.

Background

Any digital filter can be represented in its state-space format as follows: \begin{align} x(t+1) &= F x(t) + G u(t) \notag \\ y(t) &= H x(t) + D u(t) \notag \end{align}

This filter can be expressed as \begin{align} \mathcal{Y} = \mathcal{H} U + \mathcal{O}x_0 \end{align} where $\mathcal{Y} = [y_0, y_1,\ldots,y_{N-1}]^T$ is the vector of outputs, $U = [u_0,u_1,\ldots,u_{N-1}]^T$ is the vector of inputs, $x_0$ is the initial state, \begin{align} \mathcal{H} = \begin{bmatrix} D & 0 & \ldots & 0 \\ HG & D & 0 & \vdots \\ \vdots & & \ddots & \\ HF^{N-2}G & \ldots & HG & D \end{bmatrix}, \end{align} and \begin{align} \mathcal{O} = \begin{bmatrix} H \\ HF \\ \vdots \\ HF^{N-1} \end{bmatrix}. \end{align}

Problem Statement

In forward-backward filtering proposed in [1], the author mentions that the forward and backward filters are different (generally speaking). However, if we look at the implementation of the forward-backward filter derivation, the same filter is used, i.e., only the time-reversal of the input and outputs are performed but the filter transfer function does not change. What should be an efficient representation of the state-space filter while applying the forward filter and backward filter, i.e., determine $(F_f, G_f, H_f, D_f)$ and $(F_b, G_b, H_b, D_b)$.

Reference

[1] F. Gustafsson, "Determining the initial states in forward-backward filtering," in IEEE Transactions on Signal Processing, vol. 44, no. 4, pp. 988-992, Apr 1996. http://www.diva-portal.org/smash/get/diva2:315708/FULLTEXT02

Answer 1

The backwards filter, as the name suggests filters the input signal with the same filter $(F,G,H,D)$ backwards in time

$$ \begin{align} x(t-1) &= F\,x(t) + G\,u(t) \\ y(t) &= H\,x(t) + D\,u(t) \end{align} \tag{1} $$

But in order to define a filter using a state space model, the update law does have to go forwards in time, which can be done by solving the first equation for $x(t)$

$$ x(t) = F^{-1}\,\left(x(t-1) - G\,u(t)\right) \tag{2} $$

however $x(t)$ can only a function of $x(t-1)$ and $u(t-1)$. This can be solved by using a new definition for the state $z(t) = x(t) + F^{-1}\,G\,u(t)$, so

$$ z(t) = x(t) + F^{-1}\,G\,u(t) = F^{-1}\,x(t-1) = F^{-1}\,\left(z(t-1) - F^{-1}\,G\,u(t-1)\right) \tag{3} $$

which can be rewritten in the following more standard form

$$ z(t+1) = F^{-1}\,z(t) - F^{-2}\,G\,u(t). \tag{4} $$

Substituting the definition of $z(t)$ in to the expression for $y(t)$ gives

$$ y(t) = H\,x(t) + D\,u(t) = H\,\left(z(t) - F^{-1}\,G\,u(t) \right) + D\,u(t) \tag{5} $$

which can be simplified down to

$$ y(t) = H\,z(t) + \left(D - H\,F^{-1}\,G\right)\,u(t). \tag{6} $$

Combining equation $(4)$ and $(6)$ gives the following (forwards) state space model, which is equivalent to the (backwards) state space model in equation $(1)$

$$ \begin{align} z(t+1) &= F^{-1}\,z(t) - F^{-2}\,G\,u(t) \\ y(t) &= H\,z(t) + \left(D - H\,F^{-1}\,G\right)\,u(t) \end{align} \tag{7} $$

So the backwards filter based on the forwards filter $(F,G,H,D)$ can be written as $(F^{-1},-F^{-2}\,G,H,D - H\,F^{-1}\,G)$.