State-Space Representation of Forward and Backward Filters 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) &amp;= F x(t) + G u(t) \notag \\
y(t) &amp;= 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  &amp; 0 &amp; \ldots &amp; 0 \\ 
                              HG &amp; D &amp; 0 &amp; \vdots \\
                              \vdots &amp; &amp; \ddots &amp; \\
                              HF^{N-2}G &amp; \ldots &amp; HG &amp; 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
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) &amp;= F\,x(t) + G\,u(t) \\
y(t) &amp;= 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) &amp;= F^{-1}\,z(t) - F^{-2}\,G\,u(t) \\
y(t) &amp;= 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)$.