State-Space Representation – An Engineering Primer.

This is the first blog in what will probably be a series on the benefits of using state-space representation to implement discrete infinite impulse response (IIR) filters. Before I write about state-space and some of the things you can do with it, I want to give an overview of the difference equation and explain (hopefully intuitively) why they can be a poor choice for implementing a filter.

The discrete difference equation

The discrete difference equation takes an input value, a list of previous input values and a list of previous output values and combines them linearly to produce the next output value. This is codified mathematically in the following equation which can be recognised as the sum of two convolutions:

    \[ y[n] = \sum_{k=0}^{N_b} b[k] x[n-k] - \sum_{k=1}^{N_a} a[k] y[n-k] \]

Instead of the above, we could write the difference equation as two distinct difference equations chained together:

    \[ \begin{array}{@{}ll@{}} w[n] &= \sum_{k=0}^{N_b} b[k] x[n-k] \\ y[n] &= w[n] - \sum_{k=1}^{N_a} a[k] y[n-k] \end{array} \]

w[n] is the output of an all-zero filter which is used as input into an all-pole filter to produce y[n]. If we were to take the z-transforms of the two filters, we would see quite plainly that we can also swap the order of the filters as:

    \[ \begin{array}{@{}ll@{}} q[n] &= x[n] - \sum_{k=1}^{N_a} a[k] q[n-k] \\ y[n] &= \sum_{k=0}^{N_b} b[k] q[n-k] \end{array} \]

This time, q[n] is the output of an all-pole filter which is used as input into an all-zero filter to produce y[n]. If we set N=N_a=N_b, the benefit of this representation is that the two filters share exactly the same state and so an implementation needs not have delay lines for both output and input. This formulation has the name “Direct Form II” and is a popular implementation choice.

Regardless of whether we implement the difference equation directly or use Direct Form II, the implementation is easy to realise in a C program given a set of filter coefficients… so why would we want to do anything differently?

When difference equations go bad…

Let’s build a 6-th order elliptic filter and run the difference equation explicitly using the following GNU Octave code. The filter is designed with a passband edge frequency of 240 Hz at a 48000 Hz sample rate (which isn’t particularly low) with 6 dB of ripple in the passband and 80 dB of stopband attenuation.

N = 6;
IRLen = 8000;
[b, a] = ellip(N, 6, 80.0, 240.0/24000);
dp_y = zeros(IRLen, 1);
q = zeros(N, 1);
x = 1;
for i = 1:IRLen
  qv      = x - a(2:end) * q;
  dp_y(i) = b * [qv; q];
  q       = [qv; q(1:(N-1))];
  x       = 0;
end

And here is a plot of the impulse response and frequency response (which both look pretty good!):

Frequency Response of Double Precision Difference Equation - Everything looks good

Frequency Response of Double Precision Difference Equation

Oh yeah, Octave by default is using double-precision for everything. Because we like fast things, let’s specify that the states should always be stored using single precision by replacing the initialisation of q with:

q = single(zeros(N, 1));

And here is the new impulse response:

Frequency Response of Single Precision Difference Equation - it explodes into instability

Frequency Response of Single Precision Difference Equation

Bugger! Storing the states in single precision has lead to an instability in the filter. Small errors in the state are somehow being amplified by the coefficients of the difference equation and the output is growing without bound.

I’m going to try and give some intuition here: imagine that you are trying to synthesise a low-frequency discretely-sampled sinusoid that follows a slowly decaying exponential envelope. You’ve found that there is a z-transform that does exactly what you want from a z-transform table:

    \[ a^n \sin(\omega_0 n) u[n] \Leftrightarrow \frac{a z \sin(\omega_0)}{z^2 - 2az\cos(\omega_0) + a^2} \]

You decide to implement the difference equation directly using Direct Form II leading to states that are the previous outputs of the impulse response of the all-pole component of the system.

The intuition is this: the values of the states close to t=0 are all very close to zero – but after sufficient time, will eventually reach a very high-peak relative to their starting point without any external influence. Tiny errors in these small state values could grow by multiple orders of magnitude as the state evolves. This effectively happens at every zero crossing of the signal – and that is just illustrative component of the problem.

Ideally, the energy in the states of the system would only ever decrease when x[n]=0 regardless of the values of y[n]. Unfortunately, this almost never happens with a resonant system implemented using a difference equation. We should ask the question: can we find better states than previous outputs of the all-pole filter?

State-space

The Direct Form II difference equations presented earlier utilising q[n] and y[n] can be mutated into some fairly gross matrix equations:

    \[ \begin{tiny} \begin{array}{@{}ll@{}} \begin{pmatrix} q[n-N+1] \\ \vdots \\ q[n-2] \\ q[n-1] \\ q[n] \end{pmatrix} &= \begin{pmatrix} 0 & 1 & 0 & \dots & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & 1 & 0 \\ 0 & \dots & 0 & 0 & 1 \\ -a[N] & \dots & -a[3] & -a[2] & -a[1] \end{pmatrix} \begin{pmatrix} q[n-N] \\ \vdots \\ q[n-3] \\ q[n-2] \\ q[n-1] \end{pmatrix} + \begin{pmatrix} 0 \\ \vdots \\ \vdots \\ 0 \\ 1 \end{pmatrix} x[n] \\ y[n] &= \begin{pmatrix} b[N]-a[N]b[0] & \dots & b[2]-a[2]b[0] & b[1]-a[1]b[0] \end{pmatrix} \begin{pmatrix} q[n-N] \\ \vdots \\ q[n-2] \\ q[n-1] \end{pmatrix} + b[0] x[n] \end{array} \end{tiny} \]

Work through these and convince yourself they are equivalent. We can rewrite this entire system giving names to the matrices as:

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

These matrix equations are a state-space representation. The explicitly listed matrices before it are known as a “controllable canonical form” representation (for the case when b[0]=0). \mathbf{q} is the state vector and in controllable canonical form, it holds previous values of q[n] as described in the difference equations earlier.

So why would we want to do this? Who cares that we mapped a difference equation into a set of matrix equations? I’m going to list out some reasons which probably deserve their own blogs:

  1. Ability to analyse system using a vast array of matrix mathematical tools.
  2. Improve the numerical stability of the system.
  3. Preserving physical meanings of states during discretisation.
  4. Ability to design for efficient implementations.

In this blog, we’re going to look at points one and two briefly to attempt to address the problems shown earlier.

Firstly, let’s write the state equations in terms of their z-transforms so that we can build an expression that maps input to output:

    \[ \begin{array}{@{}ll@{}} z\mathbf{Q}(z) &= \mathbf{A}\mathbf{Q}(z) + \mathbf{B}\mathbf{X}(z) \\ \mathbf{Y}(z) &= \mathbf{C} \mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z) \end{array} \]

Doing some rearranging and substitution enables us to eliminate \mathbf{Q}(z):

    \[\begin{array}{@{}ll@{}} z\mathbf{Q}(z) &= \mathbf{A}\mathbf{Q}(z)+\mathbf{B}\mathbf{X}(z) \\ \left(\mathbf{I}z-\mathbf{A}\right)\mathbf{Q}(z) &= \mathbf{B}\mathbf{X}(z) \\ \mathbf{Q}(z) &= \left(\mathbf{I}z-\mathbf{A}\right)^{-1} \mathbf{B}\mathbf{X}(z) \\ \mathbf{Y}(z) &= \left( \mathbf{C}\left(\mathbf{I}z-\mathbf{A}\right)^{-1} \mathbf{B}+\mathbf{D} \right) \mathbf{X}(z) \\ &= \frac{\mathbf{C}\text{ adj}\left(\mathbf{I}z-\mathbf{A}\right) \mathbf{B}+\mathbf{D}\det\left(\mathbf{I}z-\mathbf{A}\right) }{\det\left(\mathbf{I}z-\mathbf{A}\right)} \mathbf{X}(z) \end{array}\]

The \det\left(\mathbf{I}z-\mathbf{A}\right) term in the denominator expands to the characteristic polynomial of \mathbf{A} in z and if we are trying to map a z-domain transfer function to a state-space system, it needs to equal the denominator of the transfer function. Just as the poles of a transfer function are roots of the denominator, the poles of a state-space system are the eigenvalues of the \mathbf{A} matrix. The numerator of the above expression is a matrix of polynomials in z showing how state-space is naturally able to handle multiple-input multiple-output (MIMO) systems.

Earlier, we asked can we pick better values for the states? Before getting to that, let’s answer a different question: can we pick different values for the states? The answer is: yes – and we can prove it trivially by substituting \mathbf{q}[n] with \mathbf{W}^{-1}\hat{\mathbf{q}}[n] where \mathbf{W} is some invertible square matrix:

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

This changes the coordinate system that the state vector exists in. We can pick any invertible matrix \mathbf{W} and use it to modify a system’s \mathbf{A}, \mathbf{B} and \mathbf{C} matrices keeping the system output \mathbf{y}[n] identical. The quantity \mathbf{W}\mathbf{A}\mathbf{W}^{-1} is a similarity transform on \mathbf{A} and it preserves the eigenvalues/characteristic polynomial of \mathbf{A} for all invertible \mathbf{W}.

State-space degrees of freedom

The coefficients of a difference equation are unique for a particular impulse response. The coefficients of the state-space representation for N>1 are not. In-fact, we have already demonstrated this previously by showing we can change the coordinate system of the state-vector.

Let’s go through the process of taking a 2×2 SISO system and finding the transfer function.

    \[ \begin{small} \begin{array}{@{}ll@{}} \frac{Y(z)}{X(z)} &= \frac{\mathbf{C}\text{ adj}\left(\mathbf{I}z-\mathbf{A}\right) \mathbf{B}+\mathbf{D}\det\left(\mathbf{I}z-\mathbf{A}\right) }{\det\left(\mathbf{I}z-\mathbf{A}\right)} \\ &= \frac{\begin{pmatrix}c_0 & c_1\end{pmatrix}\text{ adj}\begin{pmatrix}z-a_{00} & -a_{01} \\ -a_{10} & z-a_{11}\end{pmatrix} \begin{pmatrix}b_0 \\ b_1\end{pmatrix}+d\det\begin{pmatrix}z-a_{00} & -a_{01} \\ -a_{10} & z-a_{11}\end{pmatrix} }{\det\begin{pmatrix}z-a_{00} & -a_{01} \\ -a_{10} & z-a_{11}\end{pmatrix}} \\ &= \frac{\begin{pmatrix}c_0 & c_1\end{pmatrix}\begin{pmatrix}z-a_{11} & a_{01} \\ a_{10} & z-a_{00}\end{pmatrix} \begin{pmatrix}b_0 \\ b_1\end{pmatrix}+d\left((z-a_{00})(z-a_{11})-a_{01}a_{10}\right)}{(z-a_{00})(z-a_{11})-a_{01}a_{10}} \\ &= \frac{c_0(b_0(z-a_{11}) + b_1a_{01}) + c_1(b_0a_{10} + b_1(z-a_{00}))+d\left((z-a_{00})(z-a_{11})-a_{01}a_{10}\right)}{(z-a_{00})(z-a_{11})-a_{01}a_{10}} \\ &= \frac{dz^2+z(c_0b_0+c_1b_1-d(a_{00}+a_{11}))+c_0(b_1a_{01}-b_0a_{11})+c_1(b_0a_{10}-b_1a_{00})+d(a_{00}a_{11}-a_{01}a_{10})}{z^2-z(a_{00}+a_{11}) + a_{00}a_{11}-a_{01}a_{10}} \end{array} \end{small} \]

We can select all of these quantities to match a given second order transfer function – but how should we choose them? What makes for good coefficients? What does good mean? We could force some values to zero or one to minimize multiplies? If we want to improve numeric performance, given that the \mathbf{A} matrix is responsible for moving states forward in time, it seems reasonable to minimize the impact that errors in the state vector can have when multiplied by \mathbf{A} – this is akin to making \mathbf{A} as well conditioned as possible. This might be a good strategy for fixing our example 6th order system.

Orthogonal matrices are well-conditioned and for those deeply familiar with IIR design, you may also know that the state update for a coupled form IIR implementation is effectively determined by an orthonormal matrix that has been attenuated by some scalar. Coupled form state-updates can only be used to mimic a conjugate pole pair, but for every classical IIR design (Butterworth, Chebyshev I/II, Cauer/Elliptic) of even order (like our example 6th order problematic system), all the poles are in conjugate pairs. Start by defining \mathbf{A} as an orthogonal matrix:

    \[ \mathbf{A} = \begin{pmatrix}a_{00} & -a_{10} \\ a_{10} & a_{00} \end{pmatrix} \]

The transfer function becomes:

    \[ \begin{tiny} \frac{Y(z)}{X(z)} = \frac{dz^2+z(c_0b_0+c_1b_1-2da_{00})+c_0(-b_1a_{10}-b_0a_{00})+c_1(b_0a_{10}-b_1a_{00})+d(a_{00}^2+a_{10}^2)}{z^2-2za_{00} + a_{00}^2+a_{10}^2} \end{tiny} \]

There are two remaining free variables we can pick. Somewhat arbitrarily, we will select c_0=1 and c_1=0 to get:

    \[ \frac{Y(z)}{X(z)} = \frac{dz^2+z(b_0-2da_{00})-b_1a_{10}-b_0a_{00}+d(a_{00}^2+a_{10}^2)}{z^2-2za_{00} + a_{00}^2+a_{10}^2} \]

We now have a system of equations we can solve for the unknown matrix quantities that will create a state-space system that mimics a 2nd order resonant transfer function.

Chaining systems together

Consider we have two systems which we would like to connect together:

    \[ \begin{array}{@{}ll@{}} z\mathbf{Q}_1(z) &= \mathbf{A}_1\mathbf{Q}_1(z) + \mathbf{B}_1\mathbf{X}(z) \\ \mathbf{U}(z) &= \mathbf{C}_1 \mathbf{Q}_1(z) + \mathbf{D}_1\mathbf{X}(z) \end{array} \]

And:

    \[ \begin{array}{@{}ll@{}} z\mathbf{Q}_2(z) &= \mathbf{A}_2\mathbf{Q}_2(z) + \mathbf{B}_2\mathbf{U}(z) \\ \mathbf{Y}(z) &= \mathbf{C}_2 \mathbf{Q}_2(z) + \mathbf{D}_2\mathbf{U}(z) \end{array} \]

We can make direct substitutions to get:

    \[ \begin{array}{@{}ll@{}}  z\mathbf{Q}_1(z) &= \mathbf{A}_1\mathbf{Q}_1(z) + \mathbf{B}_1\mathbf{X}(z) \\ z\mathbf{Q}_2(z) &= \mathbf{A}_2\mathbf{Q}_2(z) + \mathbf{B}_2\mathbf{C}_1 \mathbf{Q}_1(z) + \mathbf{B}_2 \mathbf{D}_1\mathbf{X}(z) \\ \mathbf{Y}(z) &= \mathbf{C}_2 \mathbf{Q}_2(z) + \mathbf{D}_2\mathbf{C}_1 \mathbf{Q}_1(z) + \mathbf{D}_2\mathbf{D}_1\mathbf{X}(z) \end{array} \]

We can re-write the above as a set of block matrix equations:

    \[ \begin{array}{@{}ll@{}} z\begin{pmatrix}\mathbf{Q}_2(z) \\ \mathbf{Q}_1(z)\end{pmatrix} &= \begin{pmatrix} \mathbf{A}_2 & \mathbf{B}_2 \mathbf{C}_1 \\ \mathbf{0} & \mathbf{A}_1 \end{pmatrix} \begin{pmatrix}\mathbf{Q}_2(z) \\ \mathbf{Q}_1(z)\end{pmatrix} + \begin{pmatrix}\mathbf{B}_2 \mathbf{D}_1 \\ \mathbf{B}_1\end{pmatrix} \mathbf{X}(z) \\ \mathbf{Y}(z) &= \begin{pmatrix}\mathbf{C}_2 & \mathbf{D}_2 \mathbf{C}_1 \end{pmatrix}\begin{pmatrix}\mathbf{Q}_2(z) \\ \mathbf{Q}_1(z)\end{pmatrix} + \mathbf{D}_2 \mathbf{D}_1 \mathbf{X}(z) \end{array} \]

This provides a direct method for connecting systems together into one larger system. As \mathbf{A} is an upper-triangular block-diagonal matrix, the eigenvalues of the new system are the same as those of the existing systems and the \mathbf{B}_2\mathbf{C}_1 term in the connected system has no influence. Given that, can we decouple the sections completely? i.e. can we completely eliminate the \mathbf{B}_2\mathbf{C}_1 term by using a similarity transform as described previously? Let’s pose the question using block matrices:

    \[ \begin{pmatrix} \mathbf{A}_2 & \mathbf{0} \\ \mathbf{0} & \mathbf{A}_1 \end{pmatrix} \begin{pmatrix} \mathbf{W}_{00} & \mathbf{W}_{01} \\ \mathbf{W}_{10} & \mathbf{W}_{11} \end{pmatrix} = \begin{pmatrix} \mathbf{W}_{00} & \mathbf{W}_{01} \\ \mathbf{W}_{10} & \mathbf{W}_{11} \end{pmatrix} \begin{pmatrix} \mathbf{A}_2 & \mathbf{B}_2 \mathbf{C}_1 \\ \mathbf{0} & \mathbf{A}_1 \end{pmatrix} \]

    \[ \begin{array}{@{}ll@{}} \mathbf{A}_2 \mathbf{W}_{00} &= \mathbf{W}_{00} \mathbf{A}_2 \\ \mathbf{A}_1 \mathbf{W}_{10} &= \mathbf{W}_{10} \mathbf{A}_2 \\ \mathbf{A}_2 \mathbf{W}_{01} &= \mathbf{W}_{00} \mathbf{B}_2 \mathbf{C}_1 + \mathbf{W}_{01} \mathbf{A}_1 \\ \mathbf{A}_1 \mathbf{W}_{11} &= \mathbf{W}_{10} \mathbf{B}_2 \mathbf{C}_1 + \mathbf{W}_{11} \mathbf{A}_1 \end{array} \]

There is only one valid set of solutions for these equations:

    \[ \begin{array}{@{}ll@{}} \mathbf{W}_{00} &= \mathbf{I} \\ \mathbf{W}_{11} &= \mathbf{I} \\ \mathbf{W}_{10} &= \mathbf{0} \\ \mathbf{A}_2 \mathbf{W}_{01} - \mathbf{W}_{01} \mathbf{A}_1 &= \mathbf{B}_2 \mathbf{C}_1 \end{array} \]

Where the final equation is a Sylvester equation which can be solved without too much pain programatically (i.e. using a matrix inverse) if and only if \mathbf{A}_1 and \mathbf{A}_2 share no common eigenvalues. If we do this, we end up with the following system:

    \[ \begin{array}{@{}ll@{}} z\begin{pmatrix}\mathbf{Q}_2(z) \\ \mathbf{Q}_1(z)\end{pmatrix} &= \begin{pmatrix} \mathbf{A}_2 & \mathbf{0} \\ \mathbf{0} & \mathbf{A}_1 \end{pmatrix} \begin{pmatrix}\mathbf{Q}_2(z) \\ \mathbf{Q}_1(z)\end{pmatrix} + \mathbf{W}\begin{pmatrix}\mathbf{B}_2 \mathbf{D}_1 \\ \mathbf{B}_1\end{pmatrix} \mathbf{X}(z) \\ \mathbf{Y}(z) &= \begin{pmatrix}\mathbf{C}_2 & \mathbf{D}_2 \mathbf{C}_1 \end{pmatrix}\mathbf{W}^{-1}\begin{pmatrix}\mathbf{Q}_2(z) \\ \mathbf{Q}_1(z)\end{pmatrix} + \mathbf{D}_2 \mathbf{D}_1 \mathbf{X}(z) \end{array} \]

This will be better conditioned and have some other benefits in terms of efficient implementations. This block-diagonalisation we just performed is the equivalent of performing the partial-fraction decomposition on the cascaded system (which I’ve never seen a nice looking algorithmic solution for).

Bringing it home

Let’s take our problematic 6th order system, break it up into 3 second order state-space sections using the coupled-form mapping we described and evaluate the performace again:

Performance of the coupled-form state-space systems vs. the double precision difference equation

Performance of the coupled-form state-space systems vs. the double precision difference equation


There are 4 data-sets:

  1. The original double-precision difference equation results.
  2. The concatenated coupled-form systems without block-diagonalisation storing the states in single precision.
  3. The concatenated coupled-form systems without block-diagonalisation storing the states in a q15 format (15 bits of fractional precision).
  4. The concatenated coupled-form systems with block-diagonalisation storing the states in q15 format.

The performance difference between 2 and the block-diagonalised system with states in signle precision is negligible and so is not included in the plots. The performance difference between 1, 2 and 4 is practiacally identical in the passband. There is some fudging going on with the q15 plots because the precision is fixed but the range is not. I’ve tried to make things more even by:

  • 2-Norm matching the \mathbf{B} and \mathbf{C} matrices of the systems without diagonalisation.
  • 2-Norm matching each sub-section of the \mathbf{B} and \mathbf{C} matrices of the system with diagonalisation (which we can do because they are all operating independetly).
  • 2-Norm matching the \mathbf{C} matrices of the diagonalised and undiagonalised systems after the above has happened (which also requires modifying \mathbf{B}).

The performance of 1, 2 and 4 are all incredibly good. The performance of 3 shows that cascading leads to error build-up that significantly compromises the output.

These coupled-form state-space filter realisations are good down to very low frequencies. I have had success synthesising sub 10 Hz 16th order elliptic filters that work just fine with single-precision floating point. Next time, we’ll look at preservation of state meanings when we discretise an analogue system.

An attempt at an intuitive description of the QR decomposition using Householder reflectors

I really like the Householder matrix. I’ve blogged about its ability to generate an orthonormal basis containing a particular vector in a previous blog post. This blog post is a bit of a tutorial which will use the Householder matrix to perform a QR decomposition. Applying Householder reflectors to compute a QR decomposition is nothing new, but I want this blog post to attempt to provide some intuition into how the algorithm works starting from almost nothing. We’ll briefly visit inner products, matrix multiplication, the Householder matrix and then build a QR decomposition in C. I’m not going to use formal language to define these operations unless it’s necessary. I’m going to assume that you understand things like matrix multiplication and what a matrix inverse is, but not much more than that. The post will be written assuming complex numbers are being used (the C implementation will not… maybe I will add one later).

The inner product

The inner product takes two vectors and produces a scalar. You can think of a vector as being an array of length N values representing coordinates in some N-dimensional space.

    \[ \langle \mathbf{a}, \mathbf{b} \rangle = \sum^{N}_{k=1} a_k \overline{b_k} \]

Note the bar in the above indicates conjugation. If \mathbf{a} and \mathbf{b} are row vectors, we could write the above as:

    \[ \langle \mathbf{a}, \mathbf{b} \rangle = \mathbf{a} \mathbf{b}^H \]

Where \mathbf{b}^H is the transposed conjugate of \mathbf{b}. Two vectors are said to be orthogonal if they have an inner product of zero. If the inner product of a vector with itself is equal to one, the vector is said to be a unit-vector. The inner product of a vector with a unit-vector is the proportion of that vector that is pointing in the same direction as the unit-vector i.e. if \mathbf{u} is some unit vector and I define a new vector as:

    \[ \mathbf{b} = \mathbf{a} - \langle \mathbf{a}, \mathbf{u} \rangle \mathbf{u} \]

The vector \mathbf{b} is now orthogonal to \mathbf{u} i.e. \langle\mathbf{b},\mathbf{u}\rangle = 0.

There are some good geometric visualisations of these for two and three dimensional vectors… but everything continues to work in arbitrary dimensions. The above image contains two plots which both contain two perpendicular unit vectors (red and green) and an arbitrary vector (black). In both plots we can verify the red and green vectors are unit vectors by applying the inner product. We can also verify they are orthogonal by applying the inner product. If we compute the inner product of the black vector on the two unit vectors, then multiply each of the scalar results by the corresponding unit vectors used to compute them and sum the results, we should get back the black vector.

Matrix multiplications

We can define matrix multiplication using inner products of the vectors which make up the rows of the left-hand side matrix with the vectors which make up the columns of the right-hand side matrix.

    \[ \begin{tiny} \begin{pmatrix} \langle r_0, \overline{c_0} \rangle & \langle r_0, \overline{c_1} \rangle & \dots & \langle r_0, \overline{c_{N-1}} \rangle \\ \langle r_1, \overline{c_0} \rangle & \langle r_1, \overline{c_1} \rangle & \dots & \langle r_1, \overline{c_{N-1}} \rangle \\ \vdots & \vdots &  & \vdots \\ \langle r_{M-1}, \overline{c_0} \rangle & \langle r_{M-1}, \overline{c_1} \rangle & \dots & \langle r_{M-1}, \overline{c_{N-1}} \rangle \end{pmatrix} = \begin{pmatrix} \dots & r_0 & \dots \\ \dots & r_1 & \dots \\   & \vdots &  \\ \dots & r_{M-1} & \dots \end{pmatrix} \begin{pmatrix} \vdots & \vdots &  & \vdots \\ c_0 & c_1 & \dots & c_{N-1} \\ \vdots & \vdots &  & \vdots \end{pmatrix} \end{tiny} \]

We can use this to demonstrate that a matrix that is made up of rows (or columns) that are all orthogonal to each other is trivially invertible – particularly if the vectors are unit vectors, in which case the inverse of matrix \mathbf{A} is \mathbf{A}^H. To make this particularly clear:

    \[ \begin{tiny} \begin{pmatrix} \langle r_0, r_0 \rangle & \langle r_0, r_1 \rangle & \dots & \langle r_0, r_{N-1} \rangle \\ \langle r_1, r_0 \rangle & \langle r_1, r_1 \rangle & \dots & \langle r_1, r_{N-1} \rangle \\ \vdots & \vdots &  & \vdots \\ \langle r_{N-1}, r_0 \rangle & \langle r_{N-1}, r_1 \rangle & \dots & \langle r_{N-1}, r_{N-1} \rangle \end{pmatrix} = \begin{pmatrix} \dots & r_0 & \dots \\ \dots & r_1 & \dots \\   & \vdots &  \\ \dots & r_{N-1} & \dots \end{pmatrix} \begin{pmatrix} \vdots & \vdots &  & \vdots \\ \overline{r_0} & \overline{r_1} & \dots & \overline{r_{N-1}} \\ \vdots & \vdots &  & \vdots \end{pmatrix} \end{tiny} \]

If the rows are orthogonal to each other, the inner products off the main diagonal will all be zero. If the rows are unit vectors, the diagonal entries are by definition 1 and \mathbf{I} = \mathbf{A}\mathbf{A}^H.

The Householder matrix

Define a matrix \mathbf{P} as:

    \[\mathbf{P}=\mathbf{I} - 2\mathbf{v}\mathbf{v}^H\]

Where \mathbf{v} is a unit vector. Then:

    \[\begin{array}{@{}ll@{}} \mathbf{P}\mathbf{P} &= \mathbf{I}-2\left(\mathbf{v}\mathbf{v}^H+\mathbf{v}\mathbf{v}^H\right)+4\mathbf{v}\mathbf{v}^H\mathbf{v}\mathbf{v}^H \\ &= \mathbf{I}-4\mathbf{v}\mathbf{v}^H+4\mathbf{v}\left(\mathbf{v}^H\mathbf{v}\right)\mathbf{v}^H \\ & = \mathbf{I} \end{array}\]

Thereore \mathbf{P} is its own inverse (\mathbf{P}=\mathbf{P}^{-1}) and must also be unitary (\mathbf{P}=\mathbf{P}^H).

If we can work out a method to choose \mathbf{v} such that \mathbf{P} will contain a particular unit vector \mathbf{u} and we multiply \mathbf{P} by any scaled version of vector \mathbf{u}, we will get a vector which has only one non-zero entry (because all other rows of the matrix will be orthogonal to \mathbf{u}). I have described this process before in a previous blog post but will repeat it here with some more detail. Define a vector \mathbf{u} that we want the first column of \mathbf{P} to equal:

    \[ \begin{array}{@{}ll@{}} \mathbf{u} &= \left( \mathbf{I} - 2 \mathbf{v} \mathbf{v}^H \right) \mathbf{e}_1 \\ \frac{\mathbf{e}_1-\mathbf{u}}{2} &= \mathbf{v} \overline{v_1} \end{array} \]

If u_1 is real, we can break \mathbf{v} into its individual coordinates to solve for their values:

    \[ \begin{array}{@{}ll@{}} v_1 &= \sqrt{\frac{1-u_1}{2}} \\ v_2 &= \frac{-u_2}{2\sqrt{\frac{1-u_1}{2}}} \\ &\vdots \\ v_n &= \frac{-u_n}{2\sqrt{\frac{1-u_1}{2}}} \end{array} \]

Given the way the matrix is defined, the computation of the square root is unnecessary in an implementation. This is best seen if we expand out the 2 \mathbf{v} \mathbf{v}^H matrix:

    \[ \begin{array}{@{}ll@{}} 2 \mathbf{v} \mathbf{v}^H &= \begin{pmatrix} 2 v_1 \overline{v_1} & 2 v_1 \overline{v_2} & \dots & 2 v_1 \overline{v_n} \\ 2 v_2 \overline{v_1} & 2 v_2 \overline{v_2} & \dots & 2 v_2 \overline{v_n} \\ \vdots & \vdots & \ddots & \vdots \\ 2 v_n \overline{v_1} & 2 v_n \overline{v_2} & \dots & 2 v_n \overline{v_n} \end{pmatrix} \\ &= \begin{pmatrix} 1 - u_1 & -\overline{u_2} & \dots & -\overline{u_n} \\ -u_2 & \frac{u_2 \overline{u_2}}{1-u_1} & \dots & \frac{u_2 \overline{u_n}}{1-u_1} \\ \vdots & \vdots & \ddots & \vdots \\ -u_n & \frac{u_n \overline{u_2}}{1-u_1} & \dots & \frac{u_n \overline{u_n}}{1-u_1} \end{pmatrix} \end{array} \]

Define a scalar \alpha=\frac{1}{1-u_1} and a vector \mathbf{w}:

    \[\mathbf{w} = \begin{pmatrix}1 - u_1 \\ -u_2 \\ \vdots \\ -u_n \end{pmatrix} \]

Then 2 \mathbf{v} \mathbf{v}^H = \alpha \mathbf{w} \mathbf{w}^H. This is a particularly useful formulation because we rarely want to compute \mathbf{P} in its entirity and perform a matrix mutliplication. Rather, we compute:

    \[ \mathbf{P}\mathbf{A} = \left(\mathbf{I}-\alpha\mathbf{w}\mathbf{w}^H\right)\mathbf{A}=\mathbf{A}-\alpha\mathbf{w}\mathbf{w}^H\mathbf{A} \]

Which is substantially more efficient. A final necessary change is necessary: if u_1 is close to one (which implies all other coefficients of \mathbf{u} are approaching zero), \alpha will begin to approach infinity. There is a relatively simple change here which is to recognise that \mathbf{u} and -\mathbf{u} are parallel but pointing in different directions (their inner-product is -1). If we don’t care about the direction of the vector, we can change its sign to be the most numerically suitiable for the selection of \alpha i.e. if u_1 < 0 use \mathbf{u} otherwise use -\mathbf{u}.

This is actually very useful in many applications, one of which is the QR decomposition.

The QR decomposition

The QR decomposition breaks down a matrix \mathbf{A} into a unitary matrix \mathbf{Q} and an upper-triangular matrix \mathbf{R}. There are a number of ways to do this, but we are going use the Householder matrix.

    \[\mathbf{A} = \mathbf{Q}\mathbf{R}\]

Let’s start by defining a series of unitary transformations applied to \mathbf{A} as:

    \[\begin{array}{@{}ll@{}} \mathbf{X}_1&=\mathbf{P}_N \mathbf{A} \\ \mathbf{X}_{2} &= \begin{pmatrix} \mathbf{I}_{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_{N-1} \end{pmatrix} \mathbf{X}_{1} \\ \mathbf{X}_{3} &= \begin{pmatrix} \mathbf{I}_{2} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_{N-2} \end{pmatrix} \mathbf{X}_{2} \\ & \vdots \\ \mathbf{X}_{N-1} &= \begin{pmatrix} \mathbf{I}_{N-2} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_{2} \end{pmatrix} \mathbf{X}_{N-2} \end{array} \]

Where \mathbf{P}_n are all n-by-n unitary matrices. Notice that the first row of \mathbf{X}_1 will be unmodified by all the transforms that follow. The first two rows of \mathbf{X}_2 will be unmodified by all the transforms that follow – and so-on.

If we can somehow force the first row of \mathbf{P}_N to be the first column of \mathbf{A} scaled to be a unit vector (while keeping \mathbf{P}_N unitary), the first column of \mathbf{X}_1 will contain only one non-zero entry. We then set about finding a way to select \mathbf{P}_{N-1} such that the second column contains only two non-zero entries… and so on. Once the process is finished, \mathbf{X}_{N-1} is upper triangular and the unitary matrix which converts \mathbf{A} into this form can be written as:

    \[ \mathbf{Q}^H = \begin{pmatrix} \mathbf{I}_{N-2} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_{2} \end{pmatrix} \dots \begin{pmatrix} \mathbf{I}_{2} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_{N-2} \end{pmatrix}  \begin{pmatrix} \mathbf{I}_{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{P}_{N-1} \end{pmatrix} \mathbf{P}_N \]

And:

    \[ \begin{array}{@{}ll@{}} \mathbf{X}_{N-1} &= \mathbf{Q}^H\mathbf{A} \\ \mathbf{Q}\mathbf{X}_{N-1}&=\mathbf{A} \end{array} \]

When we follow the above process and the \mathbf{P}_n matricies are chosen to be Householder matrices, we have performed the QR decomposition using Householder matrices.

Follows is the C code that implements the decomposition. Because the algorithm repeatedly applies transformations to \mathbf{A} to eventually arrive at \mathbf{R}, we will design an API that operates in-place on the supplied matrix. The orthogonal matrix will be built in parallel. The N-1 iterations over the i variable, correspond to the N-1 transformation steps presented earlier. I make no claims about the computational performance of this code. The numerical performance could be improved by using a different method to compute the length of the vector (first thing that happens in the column loop).

/* Perform the QR decomposition of the given square A matrix.
 * A/R is stored and written in columns i.e.
 * a[0] a[3] a[6]
 * a[1] a[4] a[7]
 * a[2] a[5] a[8]
 *
 * q is stored in rows i.e.
 * q[0] q[1] q[2]
 * q[3] q[4] q[5]
 * q[6] q[7] q[8] */
void qr(double *a_r, double *q, int N) {
	int i, j, k;
 
	/* Initialise qt to be the identity matrix. */
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			q[i*N+j] = (i == j) ? 1.0 : 0.0;
		}
	}
 
	for (i = 0; i < N-1; i++) {
		double  norm, invnorm, alpha;
		double *ccol = a_r + i*N;
 
		/* Compute length of vector (needed to perform normalisation). */
		for (norm = 0.0, j = i; j < N; j++) {
			norm += ccol[j] * ccol[j];
		}
		norm = sqrt(norm);
 
		/* Flip the sign of the normalisation coefficient to have the
		 * effect of negating the u vector. */
		if (ccol[0] > 0.0)
			norm = -norm;
		invnorm = 1.0 / norm;
 
		/* Normalise the column
		 * Convert the column in-place into the w vector
		 * Compute alpha */
		ccol[i] = (1.0 - ccol[i] * invnorm);
		for (j = i+1; j < N; j++) {
			ccol[j] = -ccol[j] * invnorm;
		}
		alpha = 1.0 / ccol[i];
 
		/* Update the orthogonal matrix */
		for (j = 0; j < N; j++) {
			double acc = 0.0;
			for (k = i; k < N; k++) {
				acc += ccol[k] * q[k+j*N];
			}
			acc *= alpha;
			for (k = i; k < N; k++) {
				q[k+j*N] -= ccol[k] * acc;
			}
		}
 
		/* Update the upper triangular matrix */
		for (j = i+1; j < N; j++) {
			double acc = 0.0;
			for (k = i; k < N; k++) {
				acc += ccol[k] * a_r[k+j*N];
			}
			acc *= alpha;
			for (k = i; k < N; k++) {
				a_r[k+j*N] -= ccol[k] * acc;
			}
		}
 
		/* Overwrite the w vector with the column data. */
		ccol[i] = norm;
		for (j = i+1; j < N; j++) {
			ccol[j] = 0.0;
		}
	}
}

The algorithm follows almost everthing described on this page – albiet with some minor obfuscations (the \mathbf{w} vector is stored temporarily in the output matrix to avoid needing any extra storage). It’s worth mentioning that the first step of initialising the \mathbf{Q} matrix to \mathbf{I} isn’t necessary. The initialisation could be performed explicitly in the first iteration of the loop over i – but that introduces more conditional code.

Circular Convolution, Discrete Fourier Transforms and Toeplitz Matrix Multiplication

Circular Convolution

It’s fairly well known that multiplying the results of two discrete Fourier transforms and taking the inverse results in a circular convolution. This is fairly straight-forward to prove given the definitions of the DFT and its inverse.

Given the DFT:

    \[ X \left[ k \right] = \sum_{n=0}^{N-1} e^{-\frac{j 2 \pi n k}{N}} x \left[ n \right] \]

And the inverse DFT:

    \[ x \left[ n \right] = \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n k}{N}} X \left[ k \right] \]

We can derive convolution as:

    \[ \begin{array}{@{}ll@{}} y \left[ n \right] &= \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n k}{N}} \left( \sum_{m=0}^{N-1} e^{-\frac{j 2 \pi m k}{N}} a \left[ m \right] \right) \left( \sum_{l=0}^{N-1} e^{-\frac{j 2 \pi l k}{N}} b \left[ l \right] \right) \\ &= \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m \right] b \left[ l \right] \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi \left( n - l - m \right) k}{N}} \\ &= \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m N + \left( n - l \right) \mathrm{ mod } N \right] b \left[ l \right] \end{array} \]

The last step in the above recognises that the summation over k is only non-zero for certain values of l, m and n and we make a variable swap of m to attain the result. We can write the above in matrix form as:

    \[ \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{pmatrix} = \begin{pmatrix} a_0 & a_{N-1} & a_{N-2} & \dots & a_1 \\ a_1 & a_0 & a_{N-1} & & \\ a_2 & a_1 & a_0 & \ddots & \\ \vdots & & \ddots & \ddots & a_{N-1} \\ a_{N-1} & & & a_1 & a_0 \end{pmatrix} \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_{N-1} \end{pmatrix} \]

The matrix of a coefficients is a circulant matrix. Each row is a shifted copy of the preceeding row. Given that there exist \mathcal{O} \left( n \log n \right) algorithms for computing the DFT, we have shown that multiplying a vector by a circulant matrix has an efficient algorithm (note – this is only a computational reality for large N).

Circular Convolution with a Generalised DFT

Let’s redefine our DFT as:

    \[ X \left[ k \right] = \sum_{n=0}^{N-1} e^{-\frac{j 2 \pi n \left( k + \alpha \right) }{N}} x \left[ n \right] \]

Which has an inverse of:

    \[ x \left[ n \right] = \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n \left( k + \alpha \right)}{N}} X \left[ k \right] \]

This generalisation gives us some control over the boundary conditions of the DFT and hence the assumed data periodicity i.e. the DFT assumes the transformed data continues forever being repeated verbatim over and over – we can change this using \alpha. Let’s derive the convolution again:

    \[ \begin{array}{@{}ll@{}} y \left[ n \right] &= \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n \left( k + \alpha \right)}{N}} \left( \sum_{m=0}^{N-1} e^{-\frac{j 2 \pi m \left( k + \alpha \right)}{N}} a \left[ m \right] \right) \left( \sum_{l=0}^{N-1} e^{-\frac{j 2 \pi l \left( k + \alpha \right)}{N}} b \left[ l \right] \right) \\ &= \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m \right] b \left[ l \right] \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi \left( n - l - m \right) \left( k + \alpha \right)}{N}} \\ \end{array} \]

Given this, we can now draw up some matrices for various values of \alpha. Some interesting (and perhaps useful) values are \frac{1}{4}, \frac{3}{4} and \frac{1}{2}. We will restrict our attention to \frac{1}{2} which has a slightly different form once we obliterate the exponential:

    \[ Y[n] = \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m N + \left( n - l \right) \mathrm{ mod } 2 N \right] b \left[ l \right] \left( -1 \right)^m \]

We find that the matrix operation for this looks like:

    \[ \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{pmatrix} = \begin{pmatrix} a_0 & -a_{N-1} & -a_{N-2} & \dots & -a_1 \\ a_1 & a_0 & -a_{N-1} & & \\ a_2 & a_1 & a_0 & \ddots & \\ \vdots & & \ddots & \ddots & -a_{N-1} \\ a_{N-1} & & & a_1 & a_0 \end{pmatrix} \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_{N-1} \end{pmatrix} \]

This matrix is no-longer circulant; all entries to the right of the main diagonal have been negated. This convolution might not have practical value by itself, but the symmetry suggests that it might have value when combined with another.

Toeplitz Matrix Vector Multiplication

A Toeplitz matrix has the form:

    \[ \begin{pmatrix} t_0 & t_{-1} & t_{-2} & \dots & t_{-(N-1)} \\ t_1 & t_0 & t_{-1} & & \\ t_2 & t_1 & t_0 & \ddots & \\ \vdots & & \ddots & \ddots & t_{-1} \\ t_{N-1} & & & t_1 & t_0 \end{pmatrix} \]

There are efficient algorithms for performing Toeplitz matrix by vector multiplication that use a circular convolution algorithm. These algorithms end up throwing away much of the computed result (this can be seen in the previous link in the multiplication by a zero matrix). We can avoid this by using the symmetry defined in the previous section.

If we take the sum of a regular DFT convolution of A and X and the previously defined convolution of B and X, we are effectively computing the following matrix operation:

    \[ \begin{small} \begin{pmatrix} a_0+b_0 & a_{N-1}-b_{N-1} & a_{N-2}-b_{N-2} & \dots & a_1-b_1 \\ a_1+b_1 & a_0+b_0 & a_{N-1}-b_{N-1} & & \\ a_2+b_2 & a_1+b_1 & a_0+b_0 & \ddots & \\ \vdots & & \ddots & \ddots & a_{N-1}-b_{N-1} \\ a_{N-1}+b_{N-1} & & & a_1+b_1 & a_0+b_0 \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{N-1} \end{pmatrix} \end{small} \]

We can select values for a and b to create a multiplication of any matrix with equal elements on each diagonal – such as our Toeplitz matrix – using only DFTs of the length of the input data sequence. There are some interesting optimisations that can be made when the matrix is Hermitian to eliminate some of the DFTs entirely.

On the perils of cross-fading loops in organ samples

One common strategy when looping problematic organ samples is to employ a cross-fade. This is an irreversible audio modification that gradually transitions the samples leading up to the end of a loop to equal the samples that were leading into the start. The goal is to completely eliminate any sort of impulsive glitch and hopefully also create a “good spectral match”. While the ability to eliminate glitches is possible, creating a good spectral match might not be so simple.

We can think of an organ sample as being a linear combination of two signals:

  • A voiced/correlated/predictable component (choose your word) which represents the tonal part of the sample. Strictly speaking, the tonal component is not entirely predictable – there are continuous subtle variations in the speech of a pipe… but they are typically small and we will assume predictability.
  • An unvoiced/uncorrelated/unpredictable component which represents the pipe-noise of the sample.

Both of these components are necessary for realism in a sample.

The following image shows a cross-fade of two entirely correlated signals. The top picture contains the input signals and the cross-fade transition that will be used, the bottom contains the input signals with the cross-faded signal overlaid on top. There is nothing interesting about this picture: the two input signals were the same. The cross-fade transitions between two identical signals i.e. the output signal is equal to both of the input signals.

Crossfade of correlated signal

Crossfade of correlated signal

The next image shows a cross-fade of the same correlated signal with some uniform white noise overlaid on top. What is going on in the middle of output signal? It looks like it’s been attenuated a little bit and doesn’t appear to have a uniform distribution anymore.

Crossfade of correlated signal with uniform noise

Crossfade of correlated signal with uniform noise

This final image is a cross-fade purely consisting of two uniform random signals to add some more detail.

Crossfade of uniform noise

Crossfade of uniform noise

It turns out that summing two uniformly distributed random variables yields a new random variable with a triangular distribution (this is how noise for triangular dither gets generated). During the cross-fade, the distribution of the uncorrelated components of the signal is actually changed. Not only that, but the power of the signal is reduced in the middle of the transition. Here is an audio sample of two white noise streams being cross-faded over 2 seconds:

Listen for the level drop in the middle.

If the cross-fade is too long when looping a sample, it could make the pipe wind-noise duck over the duration of the cross-fade. While the effect is subtle, I have heard it in some samples and it’s not particularly natural. I suppose the TLDR advice I can give is:

  • If you can avoid cross-fading – avoid it.
  • If you cannot avoid cross-fading – make the cross-fade as short as possible (milliseconds) to avoid an obvious level-drop of noise during the transition.