This is a tutorial on how to build a digital implementation of a 2nd-order, continuously-variable filter (i.e. one where you can change the parameters runtime) that has dynamic behaviour that mimics an analogue filter. The digital model is linear (i.e. we are not going to model any transistor/diode behavior – we’ll do that in another blog). The resulting implementation is suitable for use in a digital synthesiser or equaliser. It is not attempting to be a perfect model – but it should be quite close.

The tutorial follows this rough guide:

• Difference equations and a recap on state-space.
• Description of the circuit we will be modelling (this can be skipped – it is included for completeness).
• Expressing the desired system digitally using a partially-transformed state-space model.
• Making the filter more useful
• Summary of implementation and demo video

Difference equations and state-space recap

Most EE grads know how to get a digital implementation of an analogue filter:

1. Pick a transfer function and select values for the natural frequency () [rads/second] and the dimensionless Q-factor () – these two controls are like the “frequency” and “resonance” knobs on a synth filter.

      

2. Discretise this using your favourite method (e.g. the Bilinear Transformation) to obtain a transfer function in .
3. Implement the difference equation which has coefficients trivially taken from the previous step.

More detail on those last two steps can be found in this older blog Three step digital Butterworth design.

This process works fairly well and difference equations are relatively-inexpensive to implement. But what if we want to attach and to knobs a user could twiddle (as is the case if we were building a synth)? Can we simply update the difference equation coefficients as the system processes data?

The answer is: it depends. Assuming you are using the bilinear transform, the discretisation process will always result in coefficients that are stable i.e. if the input is zero, the output will eventually decay to zero and changing the difference equation coefficients on the fly will never cause the system to end up in some state which it cannot recover from. However, the transient that the coefficient change introduces may be inappropriate for the use case. For example: if changing the coefficients can cause a large burst of energy in the output that takes thousands of samples to decay – well, that might not be okay in a synthesiser.

In the State-Space Representation – An Engineering Primer post, we went into a decent amount of detail of the benefits of working in state-space representation. The main relevant takeaway from that post is that while the transfer function representation is unique for a particular system, the state-space representation is not. A change of coordinates can be applied to the states which does not impact the steady-state response of the filter – but this change will impact the transient response to system parameter changes!

The state vector ( using terminology from the old blog) of a system derived from a difference equation consists of previous output samples from an all-pole filter. Every sample processed through the system shifts the states of the vector (discarding one element) and creates a new value. Think about this for a while as there is some intuition to be had. Consider a low-pass, very-low-frequency and under-damped (i.e. exhibits an exponentially-decaying sinusoidal impulse response) system. When an impulse is pushed in, the magnitude of all elements in the state vector slowly rise, eventually all have similar values (near the peak), then all start falling until they all have similar values again in the trough. A small error introduced in the state when all values are the same could lead the filter onto a completely different trajectory. This was demonstrated again in the state-space primer blog with the 6th order system’s inability to run with a difference equation with single-precision floats.

So how do we pick the 2nd order system to use? For a synth, the answer is probably: whatever sounds “best” as perceived by the designer. All systems will exhibit different behavior as their coefficients change. In this blog, we will chose an analogue circuit to approximately model – but that is just that: a choice.

The analogue circuit

The circuit below is what we are going to model. It is a Voltage-Controlled Voltage Source (VCVS) filter topology. The VCVS design permits the Q-factor and natural frequency to be controlled independently. I have explicitly labeled the voltages across the two capacitors because these are the states we are going to try and preserve in our digital implementation i.e. we’re not going to use the two previous all-pole output values as with a typical difference equation. By doing this, when we vary parameters, it will be like changing component values while leaving the instantaneous voltages across those capacitors the same – this is realistic, we can vary resistors by replacing them with potentiometers.

VCVS filter topology

Explaining how to apply circuit analysis is way out of scope for this blog. But one set of equations that describes the circuit in the Laplace domain is given below:

Before we go further, we will make some substitutions in the above equations to simplify further working. Define and to get:

At this point, we can do some substitution of the equations above to find i.e. the unique transfer function of the circuit:

This looks very similar to the resonant lowpass function given previosuly with and an additional gain of tacked on.

Given that there are no common components in determining and , the analogue VCVS filter topology has the ability to provide independent control of and (although, the overall gain of the system is affected by ).

We’ve stated what we want: the states to approximate the voltages across the capacitors. We know that we cannot do that if we use a standard difference equation obtained via a transfer function. So what next?

Straight to state-space

We derived the equations for the analogue model carefully in the previous section. We ensured that the terms only existed on the left hand side and only multiplied each of the states we said we were interested in. This enables us to write out the state-space model directly as follows:

If we define the state vector , the input vector and the output vector we can write the above in the usual state variable representation:

The top equation is the state equation. The bottom is the output equation. Recall that the state-space mapping is not unique. We can change the coordinates of the vector using any invertible matrix as demonstrated in the state-space primer post. But here, we have intentionally chosen the states we want to use so that our system behaves similarly to the circuit. What we will now do is apply the bilinear transform to this system in such a way that our output equation remains the same as in the continuous system – which it must if we are actually modelling the circuit. Defining the bilinear transform as:

Where is related to the sampling rate of the system we are building and can be used to control the specific frequency in the frequency response that maps exactly to the -domain transfer function (see the wikipedia page for the bilinear transform for more information). Make the substitution to the state update equation and attempt to make the left hand side :

Note that the input mixing term includes . Normally, this term would be moved into the output equation (so that it can be eliminated) and would cause the and matrices to change indicating that we have changed the states – but we don’t want this! Instead, we simply implement that operation on the input before mixing it.

The only thing that remains is to make the system causal i.e. the system has the following form:

This depends on future samples of . We multiply the state update equation by to get:

Converting this into a state-space difference-equation we get:

The above is what we implement. and are the same as what we derived for our continuous time system ( and ). and are defined as:

Making it more useful

If we stopped here, we would have the lowpass filter the circuit represented. But who just wants a lowpass filter?

Based on how we designed the system earlier, all the properties of the circuit are modeled by the and matrices. The and matrices can be thought of as buffering values available at various points in the circuit and mixing them with various gains to create the output i.e. the and matrices can be thought to not really impact the circuit – only how the output is tapped from it. If we let and and expand the analogue transfer function, we get:

We can use , and to create any numerator of our choosing. The denominator is completely determined by the matrix, but recall the list of resonant filters at the start of this post also share the same denominator. We can use this information to build multiple outputs from the same states and input sample! We can even build some smooth control to transition between them. Let’s try obtain these quantities from something a bit more readable:

Solving for , and :

• Setting , and makes the resonant lowpass filter.
• Setting , and makes the resonant highspass filter.
• Setting , and makes the resonant bandpass/bandstop filter dependent on the value of .

We could introduce a parameter to transition between lowpass and highpass as follows:

When , the system is lowpass. When , the system is highpass. When , we have a bandstop filter with a wideband attenuation of 6 dB. That’s kinda funky, but it would be cool if we could hack in some controllable bandpass/bandstop behaviour too. We can fudge this by setting:

Now when we have a controllable bandpass/bandstop filter dependent on the value of .

The controls in the red outline control the filter.

The “Sync” checkbox controls whether should be:

• If unchecked, a constant. The slider next to it will be in Hz.
• If checked, a multiple of the main oscillator frequency.

The “LP/BP/HP” slider directly controls . gets the linear gain corresponding to the value of “BP Gain”. The value of the “Q” slider is converted to as described earlier.

Summary of implementation and demo video

The implementation has 4 input parameters:

• is the continuous mode control. It is in the range [0, 1]. At 0, the system is lowpass. At 1, the system is highpass. At 0.5, the system is bandpass/bandstop.
• is the linear bandpass/bandstop gain. When and , the system is bandstop. When and , the output is the input. When and , the system is bandboost.
• the natural frequency of the filter – this will be close to where resonance occurs.
• the derived from the Q-factor.

For every sample:

1. Smooth the control parameters , , , with 1-pole smoothers to avoid instantaneous jumps which will cause glitches.
2. Form the and matrices and compute their discretised counterparts as described (yes, this is an entire bilinear transform per sample).
3. Evaluate the state update equation to form the new states.
4. Derive and using:

      

5. Evaluate the output equation to get the next output sample.

Here’s a video of it in action:

In several previous blogs (here and here), I’ve written about the orthonormal Householder matrix and how it can constructed such that it contains a particular unit-vector. A recent tweet from Eric Arnebäck posted a link to some code that was originally put together in 1999 by Tomas Akenine-Möller and John Hughes to find an orthogonal matrix which maps one unit-vector to another. This is a slightly more general problem which can be still solved in the same way using the Householder matrix. The vectors and matrices in this post are assumed to be real but the method can be modified to work with complex numbers.

In the previously linked posts, we described the method to find the Householder matrix given a unit-vector such that . Another way of wording this is: find the orthogonal matrix that contains as the first row. Now, we will generalise this to find an orthogonal matrix given a unit-vector such that where is another unit-vector. We start by defining to be the Householder reflector and attempt to solve for :

As the quantity is a scalar, the vector must be a scaled version of . Define and make the substitution into the previous equation:

Observe that has a singularity as approaches one. This can be addressed (as has been done in the previous posts) by negating the definition of the reflector (i.e. set ) and solving for again:

This time, will be defined as and the solution for becomes:

This shifts the singularity to . If the quantity is positive, this latter solution will be more numerically robust; if it is negative, the former will be more numerically robust. We will build an algorithm that performs the flip based on the sign of which results in a discontinuity of the resulting matrix when the sign changes. This might be important to consider in some use-cases.

We can avoid computing the square root because always ends up being used squared. We define where is defined to be either or leading to being defined as either or respectively. We choose which quantity to use based on the sign of . The following C listing is not efficient but provides a reference in very little code for arbitrary dimension vectors.

/* Find the orthonormal symmetric reflection matrix which
 * maps the given unit-vector x to the given unit-vector y.
 * As the matrix is symmetric, it can be interpreted as
 * either row or column major.
 *
 * This isn't efficient. */
void make_reflector(const float *y, const float *x, unsigned N, float *mat) {
	unsigned i, j;
	float proj, sgn, beta;
	for (i = 0, proj = 0.0f; i < N; i++)
		proj += y[i] * x[i];
	sgn  = (proj >= 0.0f) ? 1.0 : -1.0;
	beta = 1.0f / (proj + sgn);
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			float w_i = x[i] + y[i] * sgn;
			float w_j = x[j] + y[j] * sgn;
			mat[i*N+j] = w_i * beta * w_j;
		}
		mat[i*N+i] -= sgn;
	}
}

Because the resulting matrix is symmetric only elements need to be computed. The w vector could also be precomputed to avoid excessive repeated computation in the inner loop. The completely unrolled and branch-less 3-dimension case can be written as:

void make_reflector(float x[3], float y[3], float mat[3][3]) {
	float proj  = x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
	float sgn   = copysignf(1.0f, proj);
	float beta  = 1.0f / (proj + sgn);
	float w0    = x[0] + y[0] * sgn;
	float w1    = x[1] + y[1] * sgn;
	float w2    = x[2] + y[2] * sgn;
	float w0a   = w0 * beta;
	float w1a   = w1 * beta;
	float w2a   = w2 * beta;
	mat[0][0]   = w0a * w0 - sgn;
	mat[0][1]   = w0a * w1;
	mat[0][2]   = w0a * w2;
	mat[1][0]   = mat[0][1];
	mat[1][1]   = w1a * w1 - sgn;
	mat[1][2]   = w1a * w2;
	mat[2][0]   = mat[0][2];
	mat[2][1]   = mat[1][2];
	mat[2][2]   = w2a * w2 - sgn;
}

This looks very similar to the listings used to reflect to . None of this is new. The Householder matrix has been used for a very long time to solve this exact problem in the context of finding orthogonal matrices for performing the QR decomposition (google “QR decomposition using reflectors” to find less-verbose descriptions of what I’ve written on this page).

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:

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

is the output of an all-zero filter which is used as input into an all-pole filter to produce . 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:

This time, is the output of an all-pole filter which is used as input into an all-zero filter to produce . If we set , 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

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

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:

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 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 regardless of the values of . 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 and can be mutated into some fairly gross matrix equations:

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

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 ). is the state vector and in controllable canonical form, it holds previous values of 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:

Doing some rearranging and substitution enables us to eliminate :

The term in the denominator expands to the characteristic polynomial of in 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 matrix. The numerator of the above expression is a matrix of polynomials in 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 with where is some invertible square matrix:

This changes the coordinate system that the state vector exists in. We can pick any invertible matrix and use it to modify a system’s , and matrices keeping the system output identical. The quantity is a similarity transform on and it preserves the eigenvalues/characteristic polynomial of for all invertible .

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 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.

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 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 – this is akin to making 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 as an orthogonal matrix:

The transfer function becomes:

There are two remaining free variables we can pick. Somewhat arbitrarily, we will select and to get:

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:

And:

We can make direct substitutions to get:

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

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

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

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 and share no common eigenvalues. If we do this, we end up with the following system:

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

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 and matrices of the systems without diagonalisation.
• 2-Norm matching each sub-section of the and matrices of the system with diagonalisation (which we can do because they are all operating independetly).
• 2-Norm matching the matrices of the diagonalised and undiagonalised systems after the above has happened (which also requires modifying ).

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.