State-space representations and why should I care?

State-space representations of filters come up frequently in the topic of control systems but do not seem to be a focus when learning-about or implementing regular filter structures. Instead, the focus seems to be on transfer functions. I think this is quite a shame because state space representations both: map conveniently to implementations and provide a much more natural way to see and exploit the mathematics behind a filter. Their representation in terms of matrices opens up a gigantic amount of mathematical tools for manipulating them. Their usefulness is demonstrated in my previous post (Arbitrary polynomial-segment signal generation) and I figured would expand on some of the details in this post… although, this is really just scratching the surface.

I did some googling to see if I could figure out why state space is not completely ubiquitous and stumbled upon this DSP related blog post which appears to ponder the same thing. The author argues that “the reason is how electrical engineers think” – which is probably true – but surely the way we think has been conditioned by our education (and/or work).

I want this page to provide a bit of an overview of state space and hopefully provide some compelling reasons to use it. Let’s see how it goes…

A quick introduction.

At the core, state space representations consist of two matrix equations:

    \[   \begin{array}{@{}ll@{}}     z\hat{q} & = \mathbf{A} \hat{q} + \mathbf{B} \hat{x} \\     \hat{y} & = \mathbf{C} \hat{q} + \mathbf{D} \hat{x}   \end{array} \]

The choice of z (as opposed to s or another variable) is irrelevant – but intentional in this post as we are going to focus on discrete systems. The \mathbf{A} and \mathbf{B} matrices define how the previous N states (stored in the column vector \hat{q}) and input combine to create the next set of states. \mathbf{C} and \mathbf{D} show how the previous states and input are mixed to create the output. The above is equivalent to the difference equation below:

    \[   \begin{array}{@{}ll@{}}     \hat{q}[n] & = \mathbf{A} \hat{q}[n-1] + \mathbf{B} \hat{x}[n] \\     \hat{y}[n] & = \mathbf{C} \hat{q}[n-1] + \mathbf{D} \hat{x}[n]   \end{array} \]

Where \hat{x}, \hat{y} and \hat{q} are column vectors at time n.

When we design and implement IIRs using transfer functions and difference equations, the states are the previous N output samples. In a state space system, we can define the states as previous output samples, but there are many ways we can assign the states. I’m going to suggest that choosing previous output samples as states is actually quite a bad choice: it’s convenient and easy to understand in an implementation, but must imply that the “true” states need to be derived from previous output samples. Consider modelling an analogue electronic system with a capacitor and an inductor: somehow our model is predicting the voltage across the capacitor and current through the inductor using previous output samples. In theory, this is no problem – but on a computer with fixed precision, it can explode into instability depending on the system. This has already been widely discussed with high-order difference equations and is one of the main reasons for using cascades of biquad sections to implement high-order filters. With state space, we can quantify just how badly a system will perform by looking at the condition number of the \mathbf{A} matrix. The condition number is a reflection of how much small changes in state can be amplified when passed through the matrix which is a metric for defining the loss of precision in a process.

Going from a transfer function to state space

While the transfer function of a system is unique, the state space representation is not. There are an infinite number of ways to choose \mathbf{A}, \mathbf{B} and \mathbf{C} that will give exactly the same output.

Choosing good ways to map a transfer function into state space is outside the scope of this post, but we can create a system which is equivalent to the 2nd order direct-form system below:

    \[   \begin{array}{@{}ll@{}}     \frac{Y(z)}{X(z)} & = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} \\     y[n] & = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]   \end{array} \]

As:

    \[   \begin{array}{@{}ll@{}}     \mathbf{A} & = \left[\begin{array}{@{}cc@{}}       -a_1 & -a_2 \\       1 & 0     \end{array}\right] \\     \mathbf{B} & = \left[\begin{array}{@{}cc@{}}       1 \\       0     \end{array}\right] \\     \mathbf{C} & = \left[\begin{array}{@{}cc@{}}       b_1 - a_1 b_0 & b_2 - a_2 b_0     \end{array}\right] \\     \mathbf{D} & = b_0   \end{array} \]

The above is implemeneted in ss_siso_init() on GitHub. As stated in the comments, it’s not a particularly good mapping.

Writing a program to execute the filter is only needs to evaluate the matrix mutliplies. Notice that for a single-input single-output system (SISO), \mathbf{B} is a column vector, \mathbf{C} is a row vector and \mathbf{D} is a scalar – so we are exclusively doing vector-matrix multiplies. If the system happened to be 4th order (i.e. the \mathbf{A} matrix is 4×4) and we were targetting x86 with SSE, we could do the whole state update operation using SSE operations.

Going from state space to a transfer function

Going back to a transfer function from state space is more complicated – although for a 2nd order system (or a block diagonal system of 2×2 and 1×1 blocks), we can work out a closed-form solution by hand. The transfer function of a state space system is given as:

    \[   \begin{array}{@{}ll@{}}     H(z) & = \mathbf{C}\left(\mathbf{I} z - \mathbf{A}\right)^{-1}\mathbf{B} + \mathbf{D} \\      & = \frac{\mathbf{C} \mbox{adj} \left( \mathbf{I} z - \mathbf{A} \right) \mathbf{B}}{\left| \mathbf{I} z - \mathbf{A} \right|}   \end{array} \]

I’m going to argue that going in this direction is almost never necessary, so we aren’t going to think much about it. If we need to figure out the frequency responce, we can do exactly what we do with transfer functions and make the substitution z=e^{j w T_0} – but now we have a matrix inverse to perform.

Changing the states

Let’s make the substitution \hat{q}=\mathbf{W} \hat{v} where \mathbf{W} is some NxN invertible matrix. To remove the \mathbf{W} from the LHS of the state update, we will left-multiply both sides by \mathbf{W}^{-1} to obtain a new system in state space form.

    \[   \begin{array}{@{}ll@{}}     z \hat{v} & = \mathbf{W}^{-1} \mathbf{A} \mathbf{W} \hat{v} + \mathbf{W}^{-1} \mathbf{B} \hat{x} \\     \hat{y} & = \mathbf{C} \mathbf{W} \hat{v} + \mathbf{D} \hat{x}   \end{array} \]

We now have a system which still produces the same output, but the values of the internal states will take on different values (provided \mathbf{W} \neq \mathbf{I}). This can be particularly useful because:

  • We mentioned before that if the matrix \mathbf{A} has a high condition number, small changes in state variables can be amplified substantially during the state update. On a computer with limited precision types, this can lead to errors which propagate to the output through the \mathbf{C} matrix. By making this transformation, we might be able to reduce the condition number of the new state update matrix.
  • Executing the vector-matrix multiply for the state update increases in complexity as the order of the system increases with N^2. If we can find some way to diagonalise or block diagonalise the state update matrix, we could reduce this to a complexity increase in N.

Concatenation of state space systems

When we have the transfer functions of two systems and want to concatenate them together (i.e. send the output of one to the input of another), we can multiply the transfer functions to get the combined system. With state space, there is a bit more process but we can do a similar thing. Denote two systems:

    \[   \begin{array}{@{}ll@{}}     z\hat{q}_0 & = \mathbf{A}_0 \hat{q}_0 + \mathbf{B}_0 \hat{x} \\     \hat{w} & = \mathbf{C}_0 \hat{q}_0 + \mathbf{D}_0 \hat{x} \\     z\hat{q}_1 & = \mathbf{A}_1 \hat{q}_1 + \mathbf{B}_1 \hat{w} \\     \hat{y} & = \mathbf{C}_1 \hat{q}_1 + \mathbf{D}_1 \hat{w}   \end{array} \]

We can join these systems together to obtain:

    \[   \begin{array}{@{}ll@{}}     z\hat{q}_0 & = \mathbf{A}_0 \hat{q}_0 + \mathbf{B}_0 \hat{x} \\     z\hat{q}_1 & = \mathbf{A}_1 \hat{q}_1 + \mathbf{B}_1 \mathbf{C}_0 \hat{q}_0 + \mathbf{B}_1 \mathbf{D}_0 \hat{x} \\     \hat{y} & = \mathbf{C}_1 \hat{q}_1 + \mathbf{D}_1 \mathbf{C}_0 \hat{q}_0 + \mathbf{D}_1 \mathbf{D}_0 \hat{x}   \end{array} \]

If we now define a new column vector q that is the concatenation of the two column vectors q_0 and q_1, we can work out the matrices for the combined system. i.e.

    \[   \begin{array}{@{}ll@{}}     \hat{q} & = \left[\begin{array}{@{}cc@{}}       \hat{q}_0 \\       \hat{q}_1     \end{array}\right] \\     \mathbf{A} & = \left[\begin{array}{@{}cc@{}}       \mathbf{A}_0 & \mathbf{0} \\       \mathbf{B}_1 \mathbf{C}_0 & \mathbf{A}_1     \end{array}\right] \\     \mathbf{B} & = \left[\begin{array}{@{}cc@{}}       \mathbf{B}_0 \\       \mathbf{B}_1 \mathbf{D}_0     \end{array}\right] \\     \mathbf{C} & = \left[\begin{array}{@{}cc@{}}       \mathbf{D}_1 \mathbf{C}_0 & \mathbf{C}_1     \end{array}\right] \\     \mathbf{D} & = \mathbf{D}_1 \mathbf{D}_0   \end{array} \]

As we mentioned before, growing the size of the \mathbf{A} matrix increases complexity in N^2 unless we can somehow figure out a way to block diagonalise the matrix. If we choose to do this when we concatenate systems together, it means that we need to find some similarity transform that satisfies:

    \[   \begin{array}{@{}ll@{}}     \left[\begin{array}{@{}cc@{}}       \mathbf{A}_0 & \mathbf{0} \\       \mathbf{B}_1 \mathbf{C}_0 & \mathbf{A}_1     \end{array}\right] \left[\begin{array}{@{}cc@{}}       \mathbf{W}_{0,0} & \mathbf{W}_{0,1} \\       \mathbf{W}_{1,0} & \mathbf{W}_{1,1}     \end{array}\right] = \left[\begin{array}{@{}cc@{}}       \mathbf{W}_{0,0} & \mathbf{W}_{0,1} \\       \mathbf{W}_{1,0} & \mathbf{W}_{1,1}     \end{array}\right] \left[\begin{array}{@{}cc@{}}       \mathbf{A}_0 & \mathbf{0} \\       0 & \mathbf{A}_1     \end{array}\right]   \end{array} \]

As long as \mathbf{A}_0 and \mathbf{A}_1 share no common eigenvalues, this problem has a unique solution. Expanding out the above into distinct matrix equations, it turns out that:

    \[   \begin{array}{@{}ll@{}}     \mathbf{W}_{0,0} & = \mathbf{I} \\     \mathbf{W}_{0,1} & = \mathbf{0} \\     \mathbf{W}_{1,1} & = \mathbf{I} \\     \mathbf{A}_1 \mathbf{W}_{1,0} - \mathbf{W}_{1,0} \mathbf{A}_0 & = - \mathbf{B}_1 \mathbf{C}_0   \end{array} \]

The final equation is known as a Sylvester equation and there are closed-form methods for solving them.

I have an implementation which performs this entire diagonalising concatenation in the ss_siso_cat_diag() function on GitHub.

This is a really cool result because it means that if we construct a high-order filter using second order state space representations and none of those individual representations have common eigenvalues, we can completely block diagonalise the system i.e. the filters can be run completely in parallel! This is very similar to a partial fraction decomposition of transfer function. Furthermore, the state space representation of the block diagonalised system shows that there is an efficient implementation for running the system on a computer with vector support where the vectors are half the length of the filter order. Denote a 8th order system with a block diagonalised state transition matrix:

    \[   \begin{array}{@{}ll@{}}     \mathbf{A} & = \left[\begin{array}{@{}cccccccc@{}}       a_{0} & b_{0} & 0 & 0 & 0 & 0 & 0 & 0 \\       c_{0} & d_{0} & 0 & 0 & 0 & 0 & 0 & 0 \\       0 & 0 & a_{1} & b_{1} & 0 & 0 & 0 & 0 \\       0 & 0 & c_{1} & d_{1} & 0 & 0 & 0 & 0 \\       0 & 0 & 0 & 0 & a_{2} & b_{2} & 0 & 0 \\       0 & 0 & 0 & 0 & c_{2} & d_{2} & 0 & 0 \\       0 & 0 & 0 & 0 & 0 & 0 & a_{3} & b_{3} \\       0 & 0 & 0 & 0 & 0 & 0 & c_{3} & d_{3}     \end{array}\right] \\     \mathbf{B} & = \left[\begin{array}{@{}c@{}}       e_{0} \\       f_{0} \\       e_{1} \\       f_{1} \\       e_{2} \\       f_{2} \\       e_{3} \\       f_{3}     \end{array}\right] \\     \mathbf{C} & = \left[\begin{array}{@{}cccccccc@{}}       g_{0} & h_{0} & g_{1} & h_{1} & g_{2} & h_{2} & g_{3} & h_{3}     \end{array}\right] \\     \mathbf{D} & = i   \end{array} \]

Each lower-case letter represents some coefficient vector which will exist in memory. With the system laid out like this, we can write some C to implement the filtering operation (this is written using my cop vector abstraction which is fairly opaque):

#include "cop/cop_vec.h"
float filter(const float *coefs, float *state, float input) {
	/* load coefficients */
	v4f a      = v4f_ld(coefs + 0);
	v4f b      = v4f_ld(coefs + 4);
	v4f c      = v4f_ld(coefs + 8);
	v4f d      = v4f_ld(coefs + 12);
	v4f e      = v4f_ld(coefs + 16);
	v4f f      = v4f_ld(coefs + 20);
	v4f g      = v4f_ld(coefs + 24);
	v4f h      = v4f_ld(coefs + 28);
	float i    = coefs[32];
	/* load existing state */
	v4f q0     = v4f_ld(state + 0);
	v4f q1     = v4f_ld(state + 4);
	/* perform state update and compute output*/
	v4f inb    = v4f_broadcast(input);
	v4f new_q0 = v4f_mul(e, inb);
	v4f new_q1 = v4f_mul(f, inb);
	v4f out    = v4f_mul(g, q0);
	new_q0     = v4f_add(new_q0, v4f_mul(a, q0));
	new_q1     = v4f_add(new_q1, v4f_mul(c, q0));
	out        = v4f_add(out, v4f_mul(h, q1));
	new_q0     = v4f_add(new_q0, v4f_mul(b, q1));
	new_q1     = v4f_add(new_q1, v4f_mul(d, q1));
	v4f_st(state + 0, new_q0);
	v4f_st(state + 4, new_q1);
	/* horizontal adddition required here but I don't support it */
	return v4f_hadd_not_implemented(out) + i * input;
}

The point of including the code is to show how we can filter a sample through an 8th order system efficiently utilising vector parallelism with 8 multiply/MAC operations. Implementing this with cascaded biquads would:

  • Probably lead to pipelining issues as each biquad is dependent on the output of the previous one.
  • Probably not be implementable using vector support.

Multiple input processing

As a final point, if we wanted to derive a way to process multiple input samples through a high-order difference equation, there is a substantial amount of tedious substitution to be done. If we stay in a state space representation, there is a single substitution to make. For example:

    \[   \begin{array}{@{}ll@{}}     \hat{q}[n]        & = \mathbf{A} \hat{q}[n-1] + \mathbf{B} \hat{x}[n] \\     \hat{q}[n+1]      & = \mathbf{A} \hat{q}[n] + \mathbf{B} \hat{x}[n+1] \\                       & = \mathbf{A}^2 \hat{q}[n-1] + \mathbf{A} \mathbf{B} \hat{x}[n] + \mathbf{B} \hat{x}[n+1] \\     \hat{q}[n+\sigma] & = \mathbf{A}^{\sigma+1}\hat{q}[n-1] + \sum_{k=0}^{\sigma}\mathbf{A}^{\sigma-k}\mathbf{B}\hat{x}[n+k]  \\     \hat{y}[n]        & = \mathbf{C} \hat{q}[n-1] + \mathbf{D} \hat{x}[n] \\     \hat{y}[n+1]      & = \mathbf{C} \hat{q}[n] + \mathbf{D} \hat{x}[n+1] \\                       & = \mathbf{C} \mathbf{A} \hat{q}[n-1] + \mathbf{C} \mathbf{B} \hat{x}[n] + \mathbf{D} \hat{x}[n+1] \\     \hat{y}[n+\sigma] & = \mathbf{C} \mathbf{A}^{\sigma} \hat{q}[n-1] + \mathbf{C} \sum_{k=0}^{\sigma-1}\mathbf{A}^{\sigma-1-k}\mathbf{B} \hat{x}[n+k] + \mathbf{D} \hat{x}[n+\sigma]   \end{array} \]

The final expressions above for \hat{q} and \hat{y} are only valid for \sigma > 0.

This may or may not be useful depending on what you’re trying to do, but in my previous post I made heavy use of it to derive simple expressions for when x[n] is a polynomial in n. It would likely also be of use in building resampling filters where many output samples are discarded and most input samples are zero.

Leave a Reply

Your email address will not be published. Required fields are marked *