Virtual Analog Synthesizer: Lowpass filter 2017-01-06

During my 2nd year of studies at university we had a course on Fourier Analysis. I thoroughly enjoyed the topic and the concepts being introduced, really a super cool area! During this time I also started really digging 80s music. Some pretty fantastic songs from that era! With these two influences in my life and an electronics course, I started to get a bit obsessed with synthesizers...

(3 people playing keytars!)

I thought it would be a fun project to work on a software synthesizer which tries to emulate subtractive synthesis (the technique most commonly used by analog synthesizers). The general idea of subtractive synthesis is to have a source sound with "rich harmonic content", that is, has more harsh sounding noise such as a square-wave or a sawtooth-wave as opposed to a pure sine-wave, thus giving us more frequencies to work with that can be "sculpted" away.

The "sculpting" is performed by using different effects and stages that alter the sound. A crucial stage is filtering. Real world analog synthesizers often feature a large knob for a lowpass filter that controls how much of the higher frequencies should be removed. In analog synthesizers this filtering is realised by an electronic filter circuit. There are many circuit designs for filters that have applications in many different areas. For audio and music, the state-variable filter and the Moog-synth style 4 pole transistor ladder seem to be used often. These allow independent control over the "resonance" or Q-factor, which can give a nice sounding effect to the filtering.

If we wish to emulate an analog synthesizer, simulating a good filter is quite important. I will try to document my attempt at digitally realising a State variable filter.

Super-simple "Leaky Integrator" Filter

A simpler, but not-so-high-quality, lowpass filter can be achieved using a leaky integrator. For each new audio sample you add a version of the previous sample that has been reduced by a factor alpha to the current sample. The current sample is in turn reduced by a factor of 1 - alpha. This new value is then stored in a memory position for use when dealing with the next sample.

A higher alpha will result in less high-frequency content.

// leaky.c

float ymem1 = 0.f;

void leaky(float *value, float alpha) {
    float out = ymemm1 * alpha + (1.f - alpha) * (*value);
    ymem1 = out;
    *value = out;
}

It works quite nicely... Below is a video of a square-wave signal being treated by the leaky integrator. The alpha factor is being modulated with a low frequency sine-wave.

State Variable Filter

The state variable filter is very popular in analog music synthesis. It has the nice property that the cutoff-frequency and the Q factor can be altered independently from another. It also has the benefit of simultaneously giving outputs for a High-pass filter, Band-pass filter and a Lowpass filter, depending on which part of the circuit you tap to get an output.

The final goal is to simulate a digital version of the state variable filter. In order to do this we must first analyse the analog circuit and find the transfer function as well as the design parameters that control the Q-factor and the cutoff-frequency. After this we can transform it into its digital counterpart.

(As we will see, analysing this circuit is not necessary for us to get the transfer function for digital purposes. If we know that we want a 2 pole lowpass filter we can immediately jump to the last expression \eqref{eq:HLP}, since it has a standard form. I realised this close to the end of my derivation... But it might still be of interest to derive it.)

Analysing the Analog SVF

From the wikipedia article we can see the electrical circuit diagram for a state variable filter. The diagram can then be analysed to get a transfer function and design parameters.

circuit diagram for a State Variable filter

U2

We start by looking at the integrator circuit around op-amp U2.

Integrator op-amp circuit U2

Using the KCL rule (Kirchoff's current rule) at node $n_1$, gives us the transfer function of this part of the circuit. The voltage at $n_1$ is designated as $V$. To deal with the op-amp, we assume it is ideal, meaning that there is no voltage difference between the $+$ and $-$ poles, and since the $-$ pole is connected to ground we can set $V = 0$.

\begin{equation} \begin{aligned} \frac{V-V_{HP}}{R_f} + \frac{V-V_{BP}}{\frac{1}{j \omega C}} &= 0\newline (V-V_{BP}) R_f &= (V_{HP} - V) \frac{1}{j \omega C}\newline -V_{BP} R_f &= V_{HP} \frac{1}{j \omega C}\newline \frac{V_{BP}}{V_{HP}} &= \frac{-1}{j \omega R_f C} \end{aligned} \label{eq:VBPVHP} \end{equation}

U3

Now, looking at the integrator circuit for U3 we have pretty much the same situation.

Integrator op-amp circuit U3

Following the same procedure as for U2, we get the following transfer function

\begin{equation} \frac{V_{LP}}{V_{BP}} = \frac{-1}{j \omega R_f C} , . \label{eq:VLPVBP} \end{equation}

U1

Integrator op-amp circuit U1

Applying the KCL rule at $n_1$ and $n_2$ gives us the following two equations

\begin{equation} \begin{cases} \frac{0 - V}{R_q} + \frac{V_{BP} - V}{R_4} = 0 , , , \text{(At $n_1$)} \newline \frac{V_{in} - V}{R_g} + \frac{V_{LP} - V}{R_1} + \frac{V_{HP} - V}{R_2} = 0 , , , \text{(At $n_2$)} \end{cases} . \end{equation}

Solving this system by eliminating $V$ gives us the following

\begin{equation} \frac{1}{R_{g}} \left(- \frac{R_{q} V_{bp}}{R_{3} + R_{q}} + V_{in}\right) + \frac{1}{R_{2}} \left(- \frac{R_{q} V_{bp}}{R_{3} + R_{q}} + V_{hp}\right) + \frac{1}{R_{1}} \left(- \frac{R_{q} V_{bp}}{R_{3} + R_{q}} + V_{lp}\right) = 0 , . \label{eq:U1} \end{equation}

Getting the Lowpass Transfer Function

We want the transfer function for the lowpass filter. That is $H_{LP} = \frac{V_{LP}}{V_{in}}$. From equations \eqref{eq:VBPVHP} and \eqref{eq:VLPVBP} we can express $V_{HP}$ and $V_{BP}$ in terms of $V_{LP}$,

\begin{equation} \begin{cases} V_{BP} = -1 j \omega R_f C V_{LP} \newline V_{HP} = -1 j \omega R_f C V_{BP} = (j \omega R_f C)^2 V_{LP} \end{cases} . \label{eq:subs} \end{equation}

Now, setting $s = j \omega$ and substituting in the expressions given by \eqref{eq:subs} into the expression derived from U1 (equation \eqref{eq:U1}), and also setting $R = R_1 = R_2 = R_g$ gives us, after a lot of algebra,

\begin{equation} H_{LP} = \frac{V_{LP}}{V_{in}} = \frac{- \frac{1}{C^2 R_f^2}}{s^2 + s \underbrace{\frac{3 R_q}{C R_f (R_3 + R_q)}}_{\frac{\omega_0}{Q}} + \frac{1}{C^2 R_f^2} } = \frac{-\omega_0^2}{s^2 + s \frac{\omega_0}{Q} + \omega_0^2} , . \label{eq:HLP} \end{equation}

Writing the system transfer function in this standard form, in which the $s^2$ term is on its own in the denominator, is a standard format, from which we can identify the parameters that control the filter characteristics. In this case, the cutoff value (in radians per second $\omega_0 = f_0 2 \pi$) is given by $\omega_0 = \frac{1}{C R_f}$ and the Q-factor is given by $Q = \frac{3 R_q}{R_3 + R_q}$. Using this, we can now choose the resistor and capacitor values so that we get a filter with our requirements.

Since we are concerned with a digital filter we actually no longer need to worry about the circuit components. In fact, like mentioned in the beginning, had we known we wanted a, so called, 2-pole lowpass filter, we could have directly jumped to the last expression in \eqref{eq:HLP}.

Digital SVF

With a continuous transfer function for the analog lowpass filter $H_{LP}$ we now have to transform it and the parameters to work in the discrete time domain. We eventually want to turn this transfer function into a "difference equation" which can then be turned into an algorithm and run as a program. Key tasks here are using a "bilinear transform" and making sure that our cutoff value $\omega_0$ still behaves as we expect after this transformation. Pages 37 to 44 in "The Art of VA Filter Design" (mirror) by Vadim Zavalishin gives good information about this.

Bilinear transform

The Bilinear transform is heavily used in digital signal processing applications to convert continuous-time systems into a discrete-time systems.

I do not know much of the theory behind the Bilinear transform. But from what I gather it is a mapping ($z = e^{sT}$, where $T$ is the time between samples) that maps the frequency domain ($s$-plane) to the discrete domain ($z$-plane). The real axis of the $s$-plane (in this field $s$ is usually means as $s = j\omega$) is mapped in the $z$-plane to the unit circle $|z| = 1$.

For discretization the expression $z = e^{sT} = \frac{e^{sT/2}}{e^{-sT/2}}$ is approximated with $\approx \frac{1 + sT/2}{1 - sT/2}$. For which the inverse mapping is $s = \frac{2}{T} \frac{z-1}{z+1}$. We can get the discrete transfer function through this substitution

\begin{equation} H_d(z) = H_a \left( \frac{2}{T} \frac{z-1}{z+1} \right) , . \label{eq:bilinearsub} \end{equation}

More info here.

Pre Warping

Wikipedia and "The Art of VA Filter Design" (mirror) have more information.

The bilinear transform will distort where the frequencies that lie in the z-plane compared to the continuous case in the s-plane. Our cutoff-parameter $\omega_0$ was set in the continuous case, its value will not correspond as we want when in the discrete domain. Before the transform, in the $\text{s-plane} , -\infty{} < \omega_0 < \infty{} , \text{and in the z-plane} ,-\frac{\pi}{T} < \omega_0 < \frac{\pi}{T} , .$ Thus, the transformation of $\omega_0$ variable into the z-plane is nonlinear.

To ensure that our cutoff-parameter lands where we expect it to when doing the transform we will need to "pre warp" the parameter that we use in the continuous case. This is done as follows,

\begin{equation} \omega_{0a} = \frac{2}{T} \tan{\left( \omega \frac{T}{2} \right)} , , \label{eq:prewarp} \end{equation}

where $\omega = f_0 2 \pi$ and $f_0$ is our wanted final cutoff-frequency in Hz and $\omega_{0a}$ is our prewarped value (in rads/sec).

Getting the Difference Equation

From the last expression in \eqref{eq:HLP} we substitute in our prewarped cutoff-frequency $\omega_{0a}$ \eqref{eq:prewarp}. We then perform our bilinear transform by substituting in $s = \frac{2}{T} \frac{z-1}{z+1}$. After more algebra and cancellations we get,

$$ \frac{ Q \tan^2{\frac{T \omega}{2}} \left( 1 + 2 z^{-1} + z^{-2} \right) }{ \left(Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} \right) + z^{-1} \left( 2 Q \tan^2{\frac{T \omega}{2}} - 2 Q \right) + z^{-2} \left( Q \tan^2{\frac{T \omega}{2}} + Q - \tan{\frac{T \omega}{2}} \right) } , . $$

It is beneficial to have the first term in the denominator be $1$. Therefore we divide the numerator and denominator by $\left(Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} \right)$ and get

\begin{equation} \frac{ \frac{Q \tan^2{\frac{T \omega}{2}} }{Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} } \left( 1 + 2 z^{-1} + z^{-2}\right) }{ 1 + z^{-1} \frac{ 2 Q \tan^2{\frac{T \omega}{2}} - 2 Q }{Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} } + z^{-2} \frac{ Q \tan^2{\frac{T \omega}{2}} + Q - \tan{\frac{T \omega}{2}} }{Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} } } , . \label{eq:Hfinal} \end{equation}

A difference equation for a system can easily be derived from an expression which has the following form,

$$ \begin{split} H(z) = \frac{Y(z)}{X(z)} &= \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{a_0 + a_1 z^{-1} + a_2 z^{-2}} \newline Y(z)(a_0 + a_1 z^{-1} + a_2 z^{-2}) &= X(z) (b_0 + b_1 z^{-1} + b_2 z^{-2}) \newline a_0 Y(z) &= b_0 X(z) + b_1 z^{-1} X(z) + b_2 z^{-2} X(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) , . \end{split} $$

When transforming into the time domain from z the $z^{-k}$ factors will result in a delay of $k$ samples. This means for example that $Y(z) = X(z) + 0.5 z^{-1} X(z)$ becomes $y[n] = x[n] + 0.5 x[n-1]$.

From \eqref{eq:Hfinal} we can now read the required coefficients!

$$ \begin{cases} a_0 = 1 \newline a_1 = \frac{ 2 Q \tan^2{\frac{T \omega}{2}} - 2 Q }{Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} } \newline a_2 = \frac{ Q \tan^2{\frac{T \omega}{2}} + Q - \tan{\frac{T \omega}{2}} }{Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} } \newline b_0 = \frac{Q \tan^2{\frac{T \omega}{2}} }{Q \tan^2{\frac{T \omega}{2}} + Q + \tan{\frac{T \omega}{2}} } \newline b_1 = 2 b_0 \newline b_2 = b_0 \end{cases} $$

Testing in Python

Using python, Jupyter and the packages numpy, scipy and sympy I created a notebook that shows how the frequency-response changes according to the values of the cutoff and Q-parameters. This is shown for both the analog case and the discrete case.

Notebook is available as .html or as .ipnyb.

Putting it Into Code

We now know how to calculate the coefficients! Only thing left to do is to put it into code.

// lowpass.c

#define PI 3.141592653589793
static float Fs = 22050; //Sample rate T = 1/Fs
static float Q = 3; //ideal lp Q = 1/sqrt(2) = 0.70710678118

float xmem1 = 0;
float xmem2 = 0;
float ymem1 = 0;
float ymem2 = 0;

void lowpass(float *value, float f0) {
    //calculate coefs...
    float w0 = f0*2*PI; //in rads/sec
    float t = tan(w0/(Fs*2));
    float t2 = t*t;
    float k = Q * t2 + Q + t;

    float a1 = 2*Q*(t2-1.f)/k;
    float a2 = (Q * t2 + Q - t)/k;
    float b0 = Q*t2/k;
    float b1 = 2*b0;
    float b2 = b0;

    float x = *value;
    float out = b0*x + b1*xmem1 + b2*xmem2 - a1*ymem1 - a2*ymem2;
    xmem2 = xmem1;
    xmem1 = x;
    ymem2 = ymem1;
    ymem1 = out;
    
    *value = out;
}

Below is a video of the filter being used in my little synth project. The Q-value is set to 3 while the cutoff-frequency f0 is set to vary using a low frequency sine-wave.

Remarks

Added (2017-01-09)

There is a common form of the digital SVF IIR filter that shows up around the web. The so called Hal Chamberlin State Variable Filter.

This was first introduced by Hal Chamberlin in his book "Musical Applications of Microprocessors". This filter has the following transfer function

$$ H(z) = \frac{k_f^2 , z^{-1}}{1-(2-k_f(k_f+k_q))z^{-1} + (1-k_f k_q)z^{-2}} , , $$

where $k_f = 2 \sin{(\pi f_c T)}$ and $k_q = \frac{1}{Q}$.

This transfer function and its coefficients look a bit different from what we arrived at above. This got me a bit worried, wondering if I did something wrong. Turns out this version of the filter is not derived using the bilinear transform and uses certain approximations.

See my question at http://dsp.stackexchange.com/ for more information.

Also, made a Jupyter notebook comparing the frequency response given by

analog biquad lowpass filter
IIR with coeffecients from Robert Bristow-Johnson's EQ Cookbook
IIR with coeffecients from Hal Chamberlin's form (above)
IIR with coeffecients that we arrived at in this blog post.

available as .html or as .ipnyb.