I was a little bit premature in hoping to get into code the other day. One last dabble with some equations and we can march steadily along into some interesting code. Promise!
We've now taken a skim look at the general procedure of the BLT and an overview of the Butterworth low-pass filter design. In doing so we've used some good analytical maths and seen how the workings fall out. Very nice for blogging with the MathJax equations, but not quite so easy to map into a general purpose algorithm for code. So, let's have a look at that. It'll also explain one of the reasons for bashing away at the pole-zero representation.
Taking our analog representation in the s-plane, we can write it in terms of a product of complex roots of the polynomials of the numerator (top) and denominator (bottom), which we call the poles \( \rho\) on the bottom and zeros \( \sigma \) on the top.
$$ H_a(s) = K_a \frac {(s-\sigma_1)(s-\sigma_2)....(s-\sigma_m)}{(s-\rho_1)(s-\rho_2)....(s-\rho_m)}= K_a \frac{\prod^M_{m=1} {(s-\sigma_m)}}{\prod^N_{n=1} {(s-\rho_n)}} $$
Now, we know how to construct a filter design with poles. What if we could do a general purpose transformation of this into a similar representation in the z-plane for the digital implementation without needing to do lots of algebra. Let's have a go!
$$ H_d(z) = H_a(s) |_{s \rightarrow \frac {2}{T}\frac {(1-z^{-1})}{(1+z^{-1})}} = K_a \frac{\prod^M_{m=1} {\left(\frac {2}{T}\frac {(1-z^{-1})}{(1+z^{-1})} -\sigma_m\right)}}{\prod^N_{n=1} {\left(\frac {2}{T}\frac {(1-z^{-1})}{(1+z^{-1})} -\rho_n\right)}} $$
If we now multiply out the \( \frac {1}{T} \frac {1}{(1+z^{-1})} \) terms:
$$ H_d(z) = K_a \frac {T^N}{T^M} \frac {(1+z^{-1})^M}{(1+z^{-1})^N} \frac{\prod^M_{m=1} {\left( 2(1-z^{-1})-\sigma_m{(1+z^{-1})} \right)}}{\prod^N_{n=1}
{\left( 2(1-z^{-1}) -\rho_n{(1+z^{-1})}\right)}} $$
Simplifying (a bit)
$$ H_d(z) = K_a T^{N-M}(1+z^{-1})^{M-N} \frac{\prod^M_{m=1} {\left( 2(1-z^{-1})-\sigma_m{(1+z^{-1})} \right)}}{\prod^N_{n=1}
{\left( 2(1-z^{-1}) -\rho_n{(1+z^{-1})}\right)}} $$
Now, just taking the zeros (numerator) for the moment:
$$ \begin{align}
\left( 2(1-z^{-1})-\sigma_m{(1+z^{-1})} \right) & = (2-2z^{-1}-\sigma_mT-Tz^{-1}) \\
& = (-z^{-1}(2+\sigma_mT)+(2-\sigma_mT)) \\
& = (2+\sigma_mT)\left(z^{-1}-\frac {(2-\sigma_mT)}{(2+\sigma_mT)}\right)
\end{align}$$
Ad following the same procedure for the bottom:
$$ H_d(z) = K_a T^{N-M} \frac {\prod^M_{m=1}{(2+\sigma_mT)^M}}{\prod^N_{n=1}{(2+\rho_nT)^N}} (1+z^{-1})^{M-N} \frac{\prod^M_{m=1} {\left(z^{-1}-\frac {(2-\sigma_mT)}{(2+\sigma_mT)}\right)}}{\prod^N_{n=1}
{\left(z^{-1}-\frac {(2-\rho_nT)}{(2+\rho_nT)}\right)}} $$
What this huge monstrosity really says is that there is a large constant at the start, followed by some additional zeros matching the difference between the number of analogue complex poles and zeros followed by a set of poles and zeros transformed into the digital domain. What is really helpful is that this can then be made into an algorithm which performs some operations on the complex roots but does not need any algebraic manipulation. Ripe for code!
We'll need one additional function to go along with this to expand out the resulting factorised digital equation into a set of coefficients, but again, that should be a routine algorithm to progressively multiply out the factors.
On next to code...
Wednesday, 23 July 2014
Tuesday, 22 July 2014
Butterworth IIR Filter
Whew! That's the groundwork in place, we're now able to move onto something a bit more interesting and hopefully get back into code at some point too.
In the previous example we looked at the general normalised second order low-pass filter function:
$$ H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1} $$
When \( Q = \frac {1}{\sqrt 2} \) this is known as the Butterworth function and gives a magnitude of \(\frac {1}{\sqrt 2}\) at the passband frequency (in the normalised function where \( \omega = 1 \) ). i.e.
$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(\frac {1}{Q^2}-2)\omega^2+1)} $$
$$ |H(j\omega)|_{\omega=1}^2 = \frac {1}{1^4 + (\sqrt 2^2 - 2)1^2 + 1}= \frac {1}{\sqrt 2} $$
The magnitude response for an analogue low-pass filter can be generalised into function of \(\frac{n}{2}\) of terms of \(\omega\) in the denominator where the Q term in the 2nd order case has been brought into a general purpose coefficient \(D_n\).
$$ |H(j\omega)|^2 = \frac {1}{D_{2n}\omega^n+...+D_2\omega^2+1} $$
Note that as this is a squared function it only contains even terms, which can be written:
$$ |H(j\omega)|^2 = \frac {1}{1+D_{2n}\omega^{2n}} $$
Graphing this function shows the Butterworth response giving the maximally flat magnitude response in the passband:
This clearly shows for the prototype filter, when the passband frequency is normalised \(\frac {\omega}{\omega_p}=1\) that the magnitude response is \( |H(j\omega)|^2 = \frac {1}{2} \).
Translating this back from \(j\omega \rightarrow s \) domains gives the following transfer function:
$$ \frac {1}{1+D_{2n}\omega^{2n}}|_{\omega \rightarrow s} \rightarrow H(s)H(-s) = \frac {1}{1+(-1)^ns^{2n}} $$
Giving 2n roots (poles) of the denominator:
$$ 1+(-1)^ns^{2n}=0 $$
$$ s^{2n} = \begin{cases}
1=e^{j2k\pi} & \text{n odd} \\
-1=e^{j(2k+1)\pi} & \text{n even} \\
\end{cases} $$
Which gives us 2n poles:
$$ s_k = e^{j[\frac {(2k+n-1)}{2n}]\pi} \qquad \text{ k=1,2,....,(2n)}$$
Each of the poles have unitary magnitude and are spread evenly around the left side of the s-plane. In all cases (except n=1) there are a set of matched conjugate pairs and in odd orders an additional pole at \(s=-j\omega\).
The poles have the following polar angle and are positioned as follows on the s-plane where the circle has unitary magnitude:
The denominator polynomial for the transfer function can then be described as:
$$ H(s) = \frac {1}{D_{B}(s)} \qquad \text{where} \qquad D_{B}(s)=\prod^n_{k=1}(s-s_k) = 1+d_1s+d_2s^2+...+d_ns^n$$
Where \(D_B(s)\) is the Butterworth polynomial which can be expressed as a set of roots (\(s_k\)) or expanded form which can be found in many places and has the following values expressed in factored or non-factored form.
The factored forms show sets of quadratics each representing the complex-conjugate pairs. The polynomial terms can be calculated by:
$$ d_k = \frac { cos \left( (k-1) \frac {\pi}{2} \right)}{sin \left( \frac {k\pi}{2n} \right)}d_{k-1} \qquad \text {k=1,2,3...,n} $$
So, with these numbers we just need to apply the s-plane to z-plane mapping, using the Bilinear Transform as described before to calculate the appropriate filter. Well, we didn't quite get to code tonight. Let's see how tomorrow evening goes.
In the previous example we looked at the general normalised second order low-pass filter function:
$$ H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1} $$
When \( Q = \frac {1}{\sqrt 2} \) this is known as the Butterworth function and gives a magnitude of \(\frac {1}{\sqrt 2}\) at the passband frequency (in the normalised function where \( \omega = 1 \) ). i.e.
$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(\frac {1}{Q^2}-2)\omega^2+1)} $$
$$ |H(j\omega)|_{\omega=1}^2 = \frac {1}{1^4 + (\sqrt 2^2 - 2)1^2 + 1}= \frac {1}{\sqrt 2} $$
The magnitude response for an analogue low-pass filter can be generalised into function of \(\frac{n}{2}\) of terms of \(\omega\) in the denominator where the Q term in the 2nd order case has been brought into a general purpose coefficient \(D_n\).
$$ |H(j\omega)|^2 = \frac {1}{D_{2n}\omega^n+...+D_2\omega^2+1} $$
Note that as this is a squared function it only contains even terms, which can be written:
$$ |H(j\omega)|^2 = \frac {1}{1+D_{2n}\omega^{2n}} $$
Graphing this function shows the Butterworth response giving the maximally flat magnitude response in the passband:
This clearly shows for the prototype filter, when the passband frequency is normalised \(\frac {\omega}{\omega_p}=1\) that the magnitude response is \( |H(j\omega)|^2 = \frac {1}{2} \).
Translating this back from \(j\omega \rightarrow s \) domains gives the following transfer function:
$$ \frac {1}{1+D_{2n}\omega^{2n}}|_{\omega \rightarrow s} \rightarrow H(s)H(-s) = \frac {1}{1+(-1)^ns^{2n}} $$
Giving 2n roots (poles) of the denominator:
$$ 1+(-1)^ns^{2n}=0 $$
$$ s^{2n} = \begin{cases}
1=e^{j2k\pi} & \text{n odd} \\
-1=e^{j(2k+1)\pi} & \text{n even} \\
\end{cases} $$
Which gives us 2n poles:
$$ s_k = e^{j[\frac {(2k+n-1)}{2n}]\pi} \qquad \text{ k=1,2,....,(2n)}$$
Each of the poles have unitary magnitude and are spread evenly around the left side of the s-plane. In all cases (except n=1) there are a set of matched conjugate pairs and in odd orders an additional pole at \(s=-j\omega\).
The poles have the following polar angle and are positioned as follows on the s-plane where the circle has unitary magnitude:
The denominator polynomial for the transfer function can then be described as:
$$ H(s) = \frac {1}{D_{B}(s)} \qquad \text{where} \qquad D_{B}(s)=\prod^n_{k=1}(s-s_k) = 1+d_1s+d_2s^2+...+d_ns^n$$
Where \(D_B(s)\) is the Butterworth polynomial which can be expressed as a set of roots (\(s_k\)) or expanded form which can be found in many places and has the following values expressed in factored or non-factored form.
The factored forms show sets of quadratics each representing the complex-conjugate pairs. The polynomial terms can be calculated by:
$$ d_k = \frac { cos \left( (k-1) \frac {\pi}{2} \right)}{sin \left( \frac {k\pi}{2n} \right)}d_{k-1} \qquad \text {k=1,2,3...,n} $$
So, with these numbers we just need to apply the s-plane to z-plane mapping, using the Bilinear Transform as described before to calculate the appropriate filter. Well, we didn't quite get to code tonight. Let's see how tomorrow evening goes.
Doing the Bilinear Transform for a 2nd order filter again
Well it was fun yesterday learning about MathJax and playing around describing the bilinear transform for a 2nd order filter. As I was looking at a couple of sites last night I noticed that a good bunch of descriptions also shown the same transformation done slightly differently so as I got pencil and paper out last night out of curiosity to see how they came to the same result I thought I'd put another post so that both can be seen:
Starting again with our 2nd-order low-pass transfer function in the s-domain:
$$H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}$$
Yesterday we used the following bilinear transform with frequency warping:
$$ s \leftarrow \frac {1}{tan \frac {\omega}{2}} \frac {1-z^{-1}}{1+z^{-1}} $$
What I'd like to do today is to use a slightly different writing of this and see how the math juggling comes out:
$$ s \leftarrow \frac {1}{K} \frac {z-1}{z+1}$$
Where \(K= tan \frac {\omega}{2}\)
This is the same transform, which can be seen by multiplying by \( \frac {z^{-1}}{z^{-1}} \). The maths just plays out a little differently.
Substituting this in to the transfer function:
$$ H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1} \leftarrow \frac {1}{\frac {1}{K^2} \frac {(z-1)^2}{(z+1)^2}+\frac {1}{KQ} \frac {(z-1)}{(z+1)}+1} $$
Multiplying top and bottom by \( K^2(z+1)^2 \) gives something slightly easier:
$$ H(z) = \frac {K^2(z+1)^2}{(z-1)^2 + \frac {K}{Q}(z-1)(z+1) + K^2(z+1)^2} $$
Expanding:
$$ H(z) = \frac {K^2(z^2 + 2z + 1)}{(z^2+2z+1)+ \frac {K}{Q}z^2 - \frac {K}{Q}+K^2(z^2+2z+1)} $$
Collecting together the terms:
$$ H(z) = \frac {K^2z^2 + 2K^2z + K^2}{z^2(1+\frac{K}{Q}+K^2) + z(-2+2K^2)+(1-\frac{K}{Q}+K^2)} $$
All very easy, which gives us our coefficients:
First step then is to get the top-half functions all looking similar, so we can take \(K^2\) and using the trig identity \(K^2 = tan^2 \frac {\omega}{2} = \frac {(1-cos \omega)}{(1+cos \omega)} \). We can use the numerator in the top and the denominator in the bottom. This makes the b coefficients:
$$ a_0 = (1+cos \omega) \left( tan^2 \frac {\omega}{2} - \frac {1}{Q} tan \frac {\omega}{2} + 1 \right) $$
Using our favourite trig identities \(K = tan \frac {\omega}{2} = \frac {sin \omega}{(1+cos \omega)} \) and \(K^2 = tan^2 \frac {\omega}{2} = \frac {(1-cos \omega)}{(1+cos \omega)} \) and substituting \( 1 = \frac {(1+cos \omega)}{(1+cos \omega)} \) gives:
$$ a_0 = (1+cos \omega) \left( tan^2 \frac {omega}{2} - \frac {1}{Q} tan \frac {\omega}{2} + 1 \right) = (1+cos \omega) \left( \frac {(1-cos \omega)}{(1+cos \omega)} - \frac {1}{Q} \frac {sin \omega}{(1+cos \omega)} + \frac {(1+cos \omega)}{(1+cos \omega)} \right) $$
And crossing out terms gives:
$$ a_0 = (1-cos \omega) - \frac {sin \omega}{Q} + (1+cos \omega) = 2 - \frac {sin \omega}{Q} $$
and similarly for \(a_2\). Now taking \(a_1\):
$$ a_1 = -2 (1+cos \omega)\left(1-tan^2 \frac {\omega}{2} \right)= -2 (1+cos \omega) \left( \frac {(1+cos \omega)}{(1+cos \omega)} - \frac {(1-cos \omega)}{(1+cos \omega)} \right) = -2\left((1+cos \omega) - (1-cos \omega)\right) = -2(2cos \omega) $$
and dividing all the terms by 2 gives:
Starting again with our 2nd-order low-pass transfer function in the s-domain:
$$H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}$$
Yesterday we used the following bilinear transform with frequency warping:
$$ s \leftarrow \frac {1}{tan \frac {\omega}{2}} \frac {1-z^{-1}}{1+z^{-1}} $$
What I'd like to do today is to use a slightly different writing of this and see how the math juggling comes out:
$$ s \leftarrow \frac {1}{K} \frac {z-1}{z+1}$$
Where \(K= tan \frac {\omega}{2}\)
This is the same transform, which can be seen by multiplying by \( \frac {z^{-1}}{z^{-1}} \). The maths just plays out a little differently.
Substituting this in to the transfer function:
$$ H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1} \leftarrow \frac {1}{\frac {1}{K^2} \frac {(z-1)^2}{(z+1)^2}+\frac {1}{KQ} \frac {(z-1)}{(z+1)}+1} $$
Multiplying top and bottom by \( K^2(z+1)^2 \) gives something slightly easier:
$$ H(z) = \frac {K^2(z+1)^2}{(z-1)^2 + \frac {K}{Q}(z-1)(z+1) + K^2(z+1)^2} $$
Expanding:
$$ H(z) = \frac {K^2(z^2 + 2z + 1)}{(z^2+2z+1)+ \frac {K}{Q}z^2 - \frac {K}{Q}+K^2(z^2+2z+1)} $$
Collecting together the terms:
$$ H(z) = \frac {K^2z^2 + 2K^2z + K^2}{z^2(1+\frac{K}{Q}+K^2) + z(-2+2K^2)+(1-\frac{K}{Q}+K^2)} $$
All very easy, which gives us our coefficients:
- \( b_2 = K^2 \)
- \( b_1 = 2K^2 \)
- \( b_0 = K^2 \)
- \( a_2 = K^2+\frac{K}{Q}+1 \)
- \( a_1 = -2(1-K^2) \)
- \( a_0 = K^2-\frac{K}{Q}+1 \)
First step then is to get the top-half functions all looking similar, so we can take \(K^2\) and using the trig identity \(K^2 = tan^2 \frac {\omega}{2} = \frac {(1-cos \omega)}{(1+cos \omega)} \). We can use the numerator in the top and the denominator in the bottom. This makes the b coefficients:
- \( b_2 = 1-cos \omega \)
- \( b_1 = 2 (1-cos \omega) \)
- \( b_0 = 1-cos \omega \)
$$ a_0 = (1+cos \omega) \left( tan^2 \frac {\omega}{2} - \frac {1}{Q} tan \frac {\omega}{2} + 1 \right) $$
Using our favourite trig identities \(K = tan \frac {\omega}{2} = \frac {sin \omega}{(1+cos \omega)} \) and \(K^2 = tan^2 \frac {\omega}{2} = \frac {(1-cos \omega)}{(1+cos \omega)} \) and substituting \( 1 = \frac {(1+cos \omega)}{(1+cos \omega)} \) gives:
$$ a_0 = (1+cos \omega) \left( tan^2 \frac {omega}{2} - \frac {1}{Q} tan \frac {\omega}{2} + 1 \right) = (1+cos \omega) \left( \frac {(1-cos \omega)}{(1+cos \omega)} - \frac {1}{Q} \frac {sin \omega}{(1+cos \omega)} + \frac {(1+cos \omega)}{(1+cos \omega)} \right) $$
And crossing out terms gives:
$$ a_0 = (1-cos \omega) - \frac {sin \omega}{Q} + (1+cos \omega) = 2 - \frac {sin \omega}{Q} $$
and similarly for \(a_2\). Now taking \(a_1\):
$$ a_1 = -2 (1+cos \omega)\left(1-tan^2 \frac {\omega}{2} \right)= -2 (1+cos \omega) \left( \frac {(1+cos \omega)}{(1+cos \omega)} - \frac {(1-cos \omega)}{(1+cos \omega)} \right) = -2\left((1+cos \omega) - (1-cos \omega)\right) = -2(2cos \omega) $$
and dividing all the terms by 2 gives:
- \( b_2 = \frac {(1-cos \omega)}{2} \)
- \( b_1 = (1-cos \omega) \)
- \( b_0 = \frac {(1-cos \omega)}{2}2 \)
- \( a_2 = (1+\frac {sin \omega}{2}) \)
- \( a_1 = -2cos \omega \)
- \( a_0 = (1-\frac {sin \omega}{2}) \)
Designing IIR Digital Filters with the Bilinear Transform
As we've seen the implementation of digital filters into code is not too hard. The complicated bit is calculating the coefficients and having suitable equations for setting those coefficients so that the filter implements the amplitude and phase characteristics as the frequency changes. Getting this right needs some understanding of the underlying theory and builds a platform from which a number of other filters can be created.
I'm going to walk through the design for a general 2nd order low pass filter and use this as the basis for understanding the different concepts and poking around a bit with the theory. This is not intended in any way to be comprehensive, there are many great posts and books out there with the full detail. However, I wanted to join up the code and theory a bit and have some of fun with the equations.
There are a number of different ways of designing digital filters. The approach I'm going to use here is looking at IIR filters using the bilinear transform method. There are other transforms for designing IIR filters and Fourier methods for FIR design. I'm going to use this approach as the maths isn't too bad and it's enough to get the kind of audio filters that are interesting to play about a bit with.
The bilinear approach starts off looking at the a 'classic' analogue design as these are well established and described, then maps this from the s-domain into the z-domain to create a digital implementation. The idea is to imitate the frequency response in some way and the approach is to find a way to mathematically imitate this in the discrete domain. (I'm very thankful for that course in discrete maths many moons ago).
Starting off then, the generalised transfer function in the s-plane describing the frequency response for a second-order low-pass filter is:
$$ H(s) = \frac {g} {s^2 + ps + g} $$
where s corresponds to the complex frequency and g is the gain to make H(0) = 1.
This is the equation for a classic damped oscillating system.
What we have here is a quadratic polynomial in the denominator (the bit of the fraction below the line). The trick is to calculate the values of g and p that make this function do what we want. First of all, let's see how this function represents functions we want to use. The frequency response is described by setting /( s = j \omega /) and then to see the amplitude response we calculate the square-root of the complex-conjugates of the transfer function. What this means is putting \( s = j \omega \) and \( s = -j \omega \) into H(s) and doing the simplifying maths.
$$ |H(\omega)| = \sqrt { H(j\omega) H(-j\omega) } $$
Which for our transfer function comes out as follows:
$$ |H(\omega)|^2 = \frac {g^2}{(\omega^4+(p^2-2g)\omega^2+g^2)} $$
\(\omega\) can now be plugged into the equation to look at the behaviour of the function over different values. In this case \(\omega\) is in the value of radians, or \(0-2\pi\). We can also tweak the parameter p to get some different responses. Which looks like this:
Before moving on we're going to simplify this and normalise the gain g=1:
$$ H(s) = \frac {1} {s^2 + ps + 1} $$
$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(p^2-2)\omega^2+1)} $$
Usually p is set to \(p=\frac {1}{Q}\) giving:
$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(\frac {1}{Q^2}-2)\omega^2+1)} $$
All good. Now the bilinear transform method of getting this function into the digital domain is to substitute \( s \leftarrow \frac {1}{K} \frac {(z-1)}{(z+1)} = \frac {1}{K} \frac {(1-z^{-1})}{(1+z^{-1})}\) using whichever from suites for the maths manipulation where \( K = tan \frac {\omega}{2} \). I've covered the maths for this juggling in a different post, so we can skip it here. What it leads to for our 2nd order filter is a bi-quadratic (e.g. quadratic function in the numerator and denominator), which represents our old friend the BiQuad digital filter:
$$ H(z) = \frac {b_2 z^{-2} + b_1 z^{-1} + b_0}{a_2 z^{-2} + a_1 z^{-1} + a_0} $$
Fantastic, just what we wanted. We now have some numbers for the coefficients for our BiQuad filter implementation, remembering that \( a_0\) is usually normalised to 1.
Now we've got the basic process there's a couple of questions to consider, how do we tune this and tweak the frequency?
In the analog case this is a whole different business, however, it's simplified a little bit in the digital filter design world as the process we've used here is to design a 'prototype' filter response normalised in gain and frequency and we're then mapping this into the digital domain using the bilinear transform using frequency warping. To map our prototype onto different frequency cut-offs we simply adjust the frequency warping term as follows:
$$ K = tan \frac {\omega}{2} \leftarrow tan \left( \pi \frac {f_c}{f_s} \right)$$
Where \(f_c\) is the corner frequency and \(f_s\) is the sampling frequency, so {f_c}{f_s} can look something like \( \frac {10,000}{44,100} \) for a 10kHz corner frequency and a 44.1kHz sample rate.
Note that the bilinear transform mapping maps the frequency in the s-domain at \( \omega=1 \leftarrow f_c \) in the z-domain for our normalised prototype filter.
I'm going to walk through the design for a general 2nd order low pass filter and use this as the basis for understanding the different concepts and poking around a bit with the theory. This is not intended in any way to be comprehensive, there are many great posts and books out there with the full detail. However, I wanted to join up the code and theory a bit and have some of fun with the equations.
There are a number of different ways of designing digital filters. The approach I'm going to use here is looking at IIR filters using the bilinear transform method. There are other transforms for designing IIR filters and Fourier methods for FIR design. I'm going to use this approach as the maths isn't too bad and it's enough to get the kind of audio filters that are interesting to play about a bit with.
The bilinear approach starts off looking at the a 'classic' analogue design as these are well established and described, then maps this from the s-domain into the z-domain to create a digital implementation. The idea is to imitate the frequency response in some way and the approach is to find a way to mathematically imitate this in the discrete domain. (I'm very thankful for that course in discrete maths many moons ago).
Starting off then, the generalised transfer function in the s-plane describing the frequency response for a second-order low-pass filter is:
$$ H(s) = \frac {g} {s^2 + ps + g} $$
where s corresponds to the complex frequency and g is the gain to make H(0) = 1.
This is the equation for a classic damped oscillating system.
What we have here is a quadratic polynomial in the denominator (the bit of the fraction below the line). The trick is to calculate the values of g and p that make this function do what we want. First of all, let's see how this function represents functions we want to use. The frequency response is described by setting /( s = j \omega /) and then to see the amplitude response we calculate the square-root of the complex-conjugates of the transfer function. What this means is putting \( s = j \omega \) and \( s = -j \omega \) into H(s) and doing the simplifying maths.
$$ |H(\omega)| = \sqrt { H(j\omega) H(-j\omega) } $$
Which for our transfer function comes out as follows:
$$ |H(\omega)|^2 = \frac {g^2}{(\omega^4+(p^2-2g)\omega^2+g^2)} $$
\(\omega\) can now be plugged into the equation to look at the behaviour of the function over different values. In this case \(\omega\) is in the value of radians, or \(0-2\pi\). We can also tweak the parameter p to get some different responses. Which looks like this:
Before moving on we're going to simplify this and normalise the gain g=1:
$$ H(s) = \frac {1} {s^2 + ps + 1} $$
$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(p^2-2)\omega^2+1)} $$
Usually p is set to \(p=\frac {1}{Q}\) giving:
$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(\frac {1}{Q^2}-2)\omega^2+1)} $$
All good. Now the bilinear transform method of getting this function into the digital domain is to substitute \( s \leftarrow \frac {1}{K} \frac {(z-1)}{(z+1)} = \frac {1}{K} \frac {(1-z^{-1})}{(1+z^{-1})}\) using whichever from suites for the maths manipulation where \( K = tan \frac {\omega}{2} \). I've covered the maths for this juggling in a different post, so we can skip it here. What it leads to for our 2nd order filter is a bi-quadratic (e.g. quadratic function in the numerator and denominator), which represents our old friend the BiQuad digital filter:
$$ H(z) = \frac {b_2 z^{-2} + b_1 z^{-1} + b_0}{a_2 z^{-2} + a_1 z^{-1} + a_0} $$
Fantastic, just what we wanted. We now have some numbers for the coefficients for our BiQuad filter implementation, remembering that \( a_0\) is usually normalised to 1.
Now we've got the basic process there's a couple of questions to consider, how do we tune this and tweak the frequency?
In the analog case this is a whole different business, however, it's simplified a little bit in the digital filter design world as the process we've used here is to design a 'prototype' filter response normalised in gain and frequency and we're then mapping this into the digital domain using the bilinear transform using frequency warping. To map our prototype onto different frequency cut-offs we simply adjust the frequency warping term as follows:
$$ K = tan \frac {\omega}{2} \leftarrow tan \left( \pi \frac {f_c}{f_s} \right)$$
Where \(f_c\) is the corner frequency and \(f_s\) is the sampling frequency, so {f_c}{f_s} can look something like \( \frac {10,000}{44,100} \) for a 10kHz corner frequency and a 44.1kHz sample rate.
Note that the bilinear transform mapping maps the frequency in the s-domain at \( \omega=1 \leftarrow f_c \) in the z-domain for our normalised prototype filter.
Monday, 21 July 2014
Doing the Bilinear Transform for a 2nd order filter
Googling you can find this explained in lots of different places to varying degrees of completeness. I wanted to step this though in small detail mainly as a practice in using the MathJax notation for myself before getting onto some of the other more interesting bits.
This post will not go into any of the explanation, I'll leave that for another (possibly more useful) one.
Firstly, let's define the Laplacian s-plane (analog) transfer function that we're looking to design (a resonant low-pass filter):
$$H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}$$
We then need to use the bilinear transformation with frequency warping:
$$ s \leftarrow \frac {1}{tan \frac {\omega}{2}} \frac {1-z^{-1}}{1+z^{-1}} $$
Substituting this into the transfer function:
$$ H(z) = \frac {1}{ \frac {1}{tan^2 \frac {\omega}{2}} \frac {(1-z^{-1})^2}{(1+z^{-1})^2} + \frac {1}{Q tan \frac {\omega}{2}} \frac {1-z^{-1}}{1+z^{-1}} + 1} $$
This can be simplified by brining \( tan^2 \frac {\omega}{2} (1+z^{-1})^2\) to the top by multiplying top and bottom by the same amount giving:
$$ H(z) = \frac { tan^2 \frac {\omega}{2} (1+z^{-1})^2}{(1+z^{-1})^2 + \frac {1}{Q} tan \frac {\omega}{2} (1-z^{-1})(1+z^{-1})+tan^2 \frac {\omega}{2}(1+z^{-1})^2} $$
We can now make the following simplifications:
$$ H(z) = \frac {(1-cos \omega)(1-2z^{-1}+z^{-2})}{(1+cos \omega) (1-2z^{-1}+z^{-2})+ \frac {(1+cos \omega)}{Q} tan \frac {\omega}{2} (1-z^{-2}) + (1+cos \omega) tan^2 \frac {\omega}{2} (1-2z^{-1}+z^{-2})} $$
Collecting together the terms for \(z^{-2}\),\(z^{-1}\),1 terms to coefficients \(a_2\),\(a_1\),\(a_0\) on the bottom and \(b_2\),\(b_1\),\(b_0\) on the top:
\( a_0 \) can then be simplified by the following two steps:
$$ a_0 = (1+cos \omega) + \frac {(1+cos \omega)}{Q} \frac {sin \omega}{1+cos \omega} + (1+cos \omega) \frac {1-cos \omega}{1+cos \omega} $$
tidying up:
$$ a_0 = (1+cos \omega) + \frac {sin \omega}{Q} + (1+cos \omega) = 2 + \frac {sin \omega}{Q} = 2 (1+ \alpha) $$
where \( \alpha = \frac {sin \omega}{2Q} \)
and similarly using the same simplifications for \( a_2 \) gives \( a_2 = 2(1- \alpha) \)
Lastly for \( a_1 \) using the trig identity for \(tan^2 \frac {\omega}{2} \) one last time:
$$ a_1 = -2(1+cos \omega) + 2(1+cos \omega) \frac {(1-cos \omega)}{(1+cos \omega)} = -2 (1+cos \omega) + 2(1-cos \omega) = 2(-2cos \omega) $$
Lining up all of the coefficients and collecting the common factors:
$$ H(z) = \frac {(1-cos \omega)}{2} \frac {b_2z^{-2}+b_1z^{-1}+b_0}{a_2z^{-2}+a_1z^{-1}+a_0} $$
With the \(b_n\) values multiplied by \(\frac {(1-cos \omega)}{2}\). Which can then be normalised to set \(a_0=1\) if required.
Which gives us the same coefficients in the RBJ functions for a low-pass filter.
This post will not go into any of the explanation, I'll leave that for another (possibly more useful) one.
Firstly, let's define the Laplacian s-plane (analog) transfer function that we're looking to design (a resonant low-pass filter):
$$H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}$$
We then need to use the bilinear transformation with frequency warping:
$$ s \leftarrow \frac {1}{tan \frac {\omega}{2}} \frac {1-z^{-1}}{1+z^{-1}} $$
Substituting this into the transfer function:
$$ H(z) = \frac {1}{ \frac {1}{tan^2 \frac {\omega}{2}} \frac {(1-z^{-1})^2}{(1+z^{-1})^2} + \frac {1}{Q tan \frac {\omega}{2}} \frac {1-z^{-1}}{1+z^{-1}} + 1} $$
This can be simplified by brining \( tan^2 \frac {\omega}{2} (1+z^{-1})^2\) to the top by multiplying top and bottom by the same amount giving:
$$ H(z) = \frac { tan^2 \frac {\omega}{2} (1+z^{-1})^2}{(1+z^{-1})^2 + \frac {1}{Q} tan \frac {\omega}{2} (1-z^{-1})(1+z^{-1})+tan^2 \frac {\omega}{2}(1+z^{-1})^2} $$
We can now make the following simplifications:
- expand \({(1-z^{-1})^2} = (1-2z^{-1}+z^{-2})\) on the top and bottom
- use the trig identity \(tan^2 \frac {\omega}{2} = \frac {(1-cos \omega)}{(1+cos \omega)}\)
- simplify \( (1-z^{-1})(1+z^{-1}) = (1-z^{-2})\)
$$ H(z) = \frac {(1-cos \omega)(1-2z^{-1}+z^{-2})}{(1+cos \omega) (1-2z^{-1}+z^{-2})+ \frac {(1+cos \omega)}{Q} tan \frac {\omega}{2} (1-z^{-2}) + (1+cos \omega) tan^2 \frac {\omega}{2} (1-2z^{-1}+z^{-2})} $$
Collecting together the terms for \(z^{-2}\),\(z^{-1}\),1 terms to coefficients \(a_2\),\(a_1\),\(a_0\) on the bottom and \(b_2\),\(b_1\),\(b_0\) on the top:
- \( b_2 = (1-cos \omega)\)
- \( b_1 = -2(1-cos \omega)\)
- \( b_0 = (1-cos \omega)\)
- \( a_2 = (1+cos \omega) - \frac {(1+cos \omega)}{Q} tan \frac {\omega}{2} + (1+cos \omega)tan^2 \frac {\omega}{2}\)
- \( a_1 = -2(1+ cos \omega) + 2(1+cos \omega) tan^2 \frac {\omega}{2}\)
- \( a_0 = (1+cos \omega) + \frac {(1+cos \omega)}{Q} tan \frac {\omega}{2} + (1+cos \omega)tan^2 \frac {\omega}{2}\)
\( a_0 \) can then be simplified by the following two steps:
- using the trig identity \(tan \frac {\omega}{2} = \frac {sin \omega}{1+cos \omega}\)
- using the trig identity \(tan^2 \frac {\omega}{2} = \frac {(1-cos \omega)}{(1+cos \omega)}\) again
$$ a_0 = (1+cos \omega) + \frac {(1+cos \omega)}{Q} \frac {sin \omega}{1+cos \omega} + (1+cos \omega) \frac {1-cos \omega}{1+cos \omega} $$
tidying up:
$$ a_0 = (1+cos \omega) + \frac {sin \omega}{Q} + (1+cos \omega) = 2 + \frac {sin \omega}{Q} = 2 (1+ \alpha) $$
where \( \alpha = \frac {sin \omega}{2Q} \)
and similarly using the same simplifications for \( a_2 \) gives \( a_2 = 2(1- \alpha) \)
Lastly for \( a_1 \) using the trig identity for \(tan^2 \frac {\omega}{2} \) one last time:
$$ a_1 = -2(1+cos \omega) + 2(1+cos \omega) \frac {(1-cos \omega)}{(1+cos \omega)} = -2 (1+cos \omega) + 2(1-cos \omega) = 2(-2cos \omega) $$
Lining up all of the coefficients and collecting the common factors:
$$ H(z) = \frac {(1-cos \omega)}{2} \frac {b_2z^{-2}+b_1z^{-1}+b_0}{a_2z^{-2}+a_1z^{-1}+a_0} $$
- \(b_2 = 1\)
- \(b_1 = -2\)
- \(b_0 = 1\)
- \(a_2 = 1 + \alpha\)
- \(a_1 = -2 cos \omega\)
- \(a_0 = 1 - \alpha\)
With the \(b_n\) values multiplied by \(\frac {(1-cos \omega)}{2}\). Which can then be normalised to set \(a_0=1\) if required.
Which gives us the same coefficients in the RBJ functions for a low-pass filter.
Adding equations into Blogger
I've been enjoying working through the implementation of filters from STK into C# and putting the RBJ filters into the same code. I had a bit of a play over the weekend to make a bunch of other usual filters which I want to blog about as I get time. As I didn't want to just drop in the code without some explanation and it's been fun revisiting the theory I wanted to blog that first.
So far I've been writing the equations in Word, then getting an image capture and posting which is pretty tedious and thought that there should be a better way so had a look at that quickly this morning. There are quite a few approaches and I haven't comprehensively looked through them but jumped at the first one that seems to do the job easily enough. I've chosen to use MathJax.
A bunch of different bloggers have explained how to get this setup which is pretty easy. This one is a nice step-by-step description. I'm a bit concerned with the comment on JavaScript variables and that I have to put this into the global layout, but rather than holding off until I've looked into all the details I'm keen to get going.
The equations look like this inline \(y = \sum_{k=0}^m {a_k}{x_k}{z}^{-k}\) and \( H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}\). And can also be put standalone which displays slightly differently:
$$ y = \sum_{k=0}^m {a_k}{x_k}{z}^{-k} $$ $$H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}$$
The text for each of these is:
y = \sum_{k=0}^m {a_k}{x_k}{z}^{-k}
and
H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}
There's a good summary tutorial on the syntax here.
So far I've been writing the equations in Word, then getting an image capture and posting which is pretty tedious and thought that there should be a better way so had a look at that quickly this morning. There are quite a few approaches and I haven't comprehensively looked through them but jumped at the first one that seems to do the job easily enough. I've chosen to use MathJax.
A bunch of different bloggers have explained how to get this setup which is pretty easy. This one is a nice step-by-step description. I'm a bit concerned with the comment on JavaScript variables and that I have to put this into the global layout, but rather than holding off until I've looked into all the details I'm keen to get going.
The equations look like this inline \(y = \sum_{k=0}^m {a_k}{x_k}{z}^{-k}\) and \( H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}\). And can also be put standalone which displays slightly differently:
$$ y = \sum_{k=0}^m {a_k}{x_k}{z}^{-k} $$ $$H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}$$
The text for each of these is:
y = \sum_{k=0}^m {a_k}{x_k}{z}^{-k}
and
H(s) = \frac {1}{s^2 + \frac {s}{Q} + 1}
There's a good summary tutorial on the syntax here.
Thursday, 17 July 2014
RBJ Filter Response in HTML5
Implementing the RBJ Filters I came across this very nice little bit of HTML5 from Aike showing the magnitude response of the LPF and HPF. If I'd seen this before I wouldn't have pasted in so many images from the excel sheet and instead put some live HTML5 graphs like this in instead.
I've modified the original to put the Bandpass and Notch responses in as well. You can see below (and play!). When I get a chance I'll put the paramters values in too as this would be a neat little design tool.
I've modified the original to put the Bandpass and Notch responses in as well. You can see below (and play!). When I get a chance I'll put the paramters values in too as this would be a neat little design tool.
Subscribe to:
Posts (Atom)