First of all remember back a couple of weeks we were playing with the STK filters? The idea is that we're going to fit into that framework, so I'm going to create a new Butterworth filter called Butter and sub-class from Iir:
using System; using System.Diagnostics; using System.Numerics; using System.Collections.Generic; // Hondrou implementation of a Butterworth filter building on the STK framework // // This approach uses the techniques from Tony Fischer at York in mkfilter // http://www-users.cs.york.ac.uk/~fisher/mkfilter/ // and builds on the Java implementation by Christian d'Heureuse // www.source-code.biz, www.inventec.ch/chdh namespace AgSynth { public class Butter : Iir { } }
First of all, it's important to give credit when building on the work of others. In this case the theoretical approach follows that described in a C implementation called mkfilter by Tony Fischer at York. More closely we're going to build some C# code flavoured by Christian d'Heureuse's Java implementation of this approach. We'll depart from both implementations by trying to put the functions into a set of more general functions that can be used to build further filter designs in the spirit of STK which was developed to understand audio synthesis and digital filtering rather than make optimised code.
Second point to note, we're using the C# System.Numerics library to bring in the Complex class for implementing some of the maths. If you haven't already you'll need to import the reference into your project (it works under Silverlight too). The really good thing here is that the 'out-of-the-box' Complex class pretty much matches Christian's Complex class. Where there are differences it's pretty easy to resolve.
So, first function, let's jump in at the deep end and make a method to generate those Complex poles for the Butterworth polynomial.
public List<Complex> getPoles(int order) { Complex[] poles = new Complex[order]; for (int i = 0; i < order; i++) { double theta = (order / 2.0 + 0.5 + i) * Math.PI / order; poles[i] = Complex.FromPolarCoordinates(1, theta); } return new List<Complex>(poles); }
Pretty easy! I'm going to continue to use List<Complex> rather than Complex[] for no other particular reason than I used List<double> for the a,b coefficients in the Filter class. It means a couple of times we'll need to use Complex[] to set arrays, but all the methods will use List<Complex>.
It's worth checking this, so using some test code below to dump out the list of Complex numbers:
private void Dump(string message, List<Complex> cs, bool polar) { String text = message; foreach (Complex c in cs) { if (polar) { text += String.Format("[r={0:0.00} p={1:0.00}] ", c.Magnitude, c.Phase); } else { text += String.Format("{0:0.00}+i{1:0.00} ", c.Real, c.Magnitude); } } Debug.WriteLine(text); } private void test() { Butter b = new Butter(); Dump("n=1 ",b.getPoles(1),false); Dump("n=2 ",b.getPoles(2),false); Dump("n=3 ",b.getPoles(3),false); Dump("n=4 ",b.getPoles(4),false); }
Gives this result:
n=1 -1.00+i1.00
n=2 -0.71+i1.00 -0.71+i1.00
n=3 -0.50+i1.00 -1.00+i1.00 -0.50+i1.00
n=4 -0.38+i1.00 -0.92+i1.00 -0.92+i1.00 -0.38+i1.00
Or in polar format:
n=1 [r=1.00 p=3.14]
n=2 [r=1.00 p=2.36] [r=1.00 p=-2.36]
n=3 [r=1.00 p=2.09] [r=1.00 p=3.14] [r=1.00 p=-2.09]
n=4 [r=1.00 p=1.96] [r=1.00 p=2.75] [r=1.00 p=-2.75] [r=1.00 p=-1.96]
Where the angle is in radians. Which matches what we would expect.
Now, let's create a new function in the Butter class to make a low-pass filter with a specified corner frequency (fcf) and order. This is the design function that we'll be working through and then building on so I'm going to do it bit-by-bit. To be consistent with the rest of the classes fcf is specified in Hz referenced to the SampleRate (which is 44100). We're going to first calculate the normalised corner frequency, check it and then use it to pre-warp the poles in preparation for a bilinear transform.
// fcf is the cutoff frequency in Hz public void lpf(double fcf, int order) { double nfcf = fcf / PcmStreamSource.SampleRate; // normalised corner frequency Debug.Assert(nfcf > 0 && nfcf <= 0.5); // greater than 0 and less than Nyquist List<Complex> poles = getPoles(order); // get warped freq double w1 = 2 * Math.PI * Math.Tan(Math.PI * nfcf) / Math.PI; List<Complex> warpedpoles = Polynomial.multiply(poles, w1); // ..... more to come
The Polynomial class is a helper class that we'll progressively add functions to. While it's tempting to make this a class that replaces List<Complex> it does not paricularly add much in aesthetic and makes a whole lot more code, so we're going to use this as a class with a bunch of static helper functions that are aimed at working with a list of complex numbers which represent a polynomial. I'm not really going to distinguish here whether they are a list of coefficients or roots.
Multiply is pretty easy, it scales each of the complex numbers:
public class Polynomial { public static List<Complex> multiply(List<Complex> a, double d) { Complex[] result = new Complex[a.Count]; for (int i = 0; i < a.Count; i++) { result[i] = a[i] * d; } return new List<Complex>(result); } }
We'll add to this class progressively as the functions are needed.
Great! So we now have a set of pre-warped poles that effectively represent our s-plane, so I'm going to create a class to represent the ComplexPlane:
public class ComplexPlane { public List<Complex> poles; public List<Complex> zeros; }
And we'll update the lpf function to create it. In the case of a lpf as we've seen, the transfer function in the s-plane does not have any zeros, so let's plug that in as well. So, from the pre-warping it now looks like this:
// get warped freq double w1 = 2 * Math.PI * Math.Tan(Math.PI * nfcf) / Math.PI; ComplexPlane splane = new ComplexPlane(); splane.poles = Polynomial.multiply(poles,w1); splane.zeros = new List<Complex>(); ComplexPlane zplane = splane.BilinearTransform();
Where we've already added the next step, which is a function in our ComplexPlane class to transform it and generate another ComplexPlane object which is the zplane. The BilinearTransform method is below, and as I promised to not do more maths in this post, please look at the previous post explaining the method and we'll maybe come back to this again if needed. For now, we get a set of poles + zeros in the zplane.
First of all, here's the function in ComplexPlane:
public ComplexPlane BilinearTransform() { ComplexPlane zplane = new ComplexPlane(); zplane.poles = Polynomial.BilinearTransform(poles); List<Complex> a = Polynomial.BilinearTransform(zeros); // extend zeros with -1 to the number of poles zplane.zeros = Polynomial.extend(a, poles.Count, new Complex(-1,0)); return zplane; }
Which requires these supporting functions in Polynomial:
public static List<Complex> BilinearTransform(List<Complex> a) { Complex[] a2 = new Complex[a.Count]; for (int i = 0; i < a.Count; i++) { a2[i] = doBilinearTransform(a[i]); } return new List<Complex>(a2); } private static Complex doBilinearTransform(Complex x) { return ((2.0 + x) / (2.0 - x)); }
And the extend function which adds a number of Complex roots to an existing root (pole or zero) representation of a factored polynomial. In this case the Bilinear Transformation adds a number of zeros:
public static List<Complex> extend(List<Complex> a, int n, Complex fill) { List<Complex> result = null; if (a.Count >= n) { result = a; } else { Complex[] a2 = new Complex[n]; for (int i = 0; i < a.Count; i++) { a2[i] = a[i]; } for (int i = a.Count; i < n; i++) { a2[i] = fill; } result = new List<Complex>(a2); } return result; }
Wow! That's a lot easier than the maths. We're now on the downhill run.
From the z-plane root representation we really want to move this into a transfer function and calculate some coefficients, so we're going to create a new class TransferFunction and be able to transform a ComplexPlane object into a TransferFunction object. The Transfer function in this case only has real coefficient values in comparison to the Complex roots in the ComplexPlane.
public class TransferFunction { public List<double> numerator; public List<double> denominator; }
And a function in the ComplexPlane object to generate a TransferFunction object:
public TransferFunction getTransferFunction() { List<Complex> numeratorComplexCoeff = Polynomial.expand(zeros); List<Complex> denominatorComplexCoeff = Polynomial.expand(poles); TransferFunction tf = new TransferFunction(); double eps = 1e-10; tf.numerator = Polynomial.toDouble(numeratorComplexCoeff, eps); tf.denominator = Polynomial.toDouble(denominatorComplexCoeff, eps); return tf; }
This uses two supporting functions in Polynomial. First expand, which takes the roots of a polynomial and expands them into a list of Complex coefficients of a polynomial.
public static List<Complex> expand(List<Complex> zeros) { List<Complex> result = null; int n = zeros.Count; if (n == 0) { result = new List<Complex>() { Complex.One }; } else { List<Complex> a = new List<Complex>() { Complex.One, Complex.Negate(zeros[0]) }; for (int i = 1; i < n; i++) { List<Complex> a2 = new List<Complex> { Complex.One, Complex.Negate(zeros[i]) }; a = multiply(a, a2); } result = a; } return result; }
The expand function itself makes use of a function of Polynomial to multiple two List<Complex> together:
public static List<Complex> multiply(List<Complex> a1, List<Complex> a2) { int n1 = a1.Count - 1; int n2 = a2.Count - 1; int n3 = n1 + n2; //List<Complex> a3 = new List<Complex>(n3+1); Complex[] a3 = new Complex[n3 + 1]; for (int i = 0; i <= n3; i++) { Complex t = Complex.Zero; int p1 = Math.Max(0, i - n2); int p2 = Math.Min(n1, i); for (int j = p1; j <= p2; j++) { Complex tmul = a1[n1 - j] * a2[n2 - i + j]; t = t + tmul; } a3[n3 - i] = t; } return new List<Complex>(a3); }
As this expansion function is pretty interesting and important in itself, it's worth checking to see that it is doing what we'd expect. So, let's take the factored representation of the poles before pre-warping and expand them as a check which should give us the Butterworth polynomial coefficients:
Dump("n=1 ", Polynomial.expand(b.getPoles(1)), true); Dump("n=2 ", Polynomial.expand(b.getPoles(2)), true); Dump("n=3 ", Polynomial.expand(b.getPoles(3)), true); Dump("n=4 ", Polynomial.expand(b.getPoles(4)), true);
This gives the following polar values:
n=1 [r=1.00 p=0.00] [r=1.00 p=0.00]
n=2 [r=1.00 p=0.00] [r=1.41 p=0.00] [r=1.00 p=0.00]
n=3 [r=1.00 p=0.00] [r=2.00 p=0.00] [r=2.00 p=0.00] [r=1.00 p=0.00]
n=4 [r=1.00 p=0.00] [r=2.61 p=0.00] [r=3.41 p=0.00] [r=2.61 p=0.00] [r=1.00 p=0.00]
Which correspond with the polynomial coefficients described previously. Great! It's worth noting, as we would hope that the angle in all cases is zero meaning that these are only real coefficients, which leads us on to the second function which converts the Complex polynomial coefficients into real values. The detail here is to compensate for a slight difficulty in the trig function implementation in the standard library where Math.Sin(Math.PI) does not quite return 0. The same is true for Excel by the way, just try it! SIN(PI()) gives the value 1.22515E-16.
// Converts complex array into real double array, verifies that all imaginary parts // are equal or below eps</code. public static List<double> toDouble(List<Complex> a, double eps) { double[] a2 = new double[a.Count]; for (int i = 0; i < a.Count; i++) { double absIm = Math.Abs(a[i].Imaginary); Debug.Assert(absIm < eps && absIm < Math.Abs(a[i].Real * eps)); a2[i] = a[i].Real; } return new List<double>(a2); }
We're very, very nearly there. Now we have our transfer function in the z domain with real coefficients we need to deal with the huge constant created in our approach to the Bilinear Transform. The easiest way to do this is to side-step it completely. The way to do this is to consider that we will want to normalise the amplitude at a particular frequency. We can then use this to scale the coefficients appropriately. So, let's add a function to TransferFunction to enable the calculation of gain at a particular frequency.
public double getGainAt(Complex w) { Debug.Assert(numerator.Count > 0); Complex sumn = new Complex(numerator[0], 0); for (int i = 1; i < numerator.Count; i++) { sumn = (sumn * w) + numerator[i]; } Debug.Assert(numerator.Count > 0); Complex sumd = new Complex(denominator[0], 0); for (int i = 1; i < denominator.Count; i++) { sumd = (sumd * w) + denominator[i]; } return (sumn / sumd).Magnitude; }
Note that as this is a z-domain transfer function in the z-plane the frequency is given as a Complex number, which is a rotating frequency vector. e.g. in the z-plane the \( \pi \) radians is the Nyquist frequency and corresponds to a vector of -1 (which is why we only consider points in the top half of the plane as frequencies in the bottom half of the plane are above the Nyquist frequency, it also explains the aliasing symmetry need to have conjugate pairs). Anyway, moving swiftly on...
And finally we can add a function to get the filter design coefficients from the transfer function, which effectively is returning a scaled version of the same transfer function.
public class Coefficients { public List<double> a; public List<double> b; }
And the function of TransferFunction:
public Coefficients getCoefficients() { double scale = denominator[0]; Coefficients coeffs = new Coefficients(); coeffs.a = new List<double>(); coeffs.b = new List<double>(); coeffs.a.Add(1.0); // a[0] for (int i = 1; i < denominator.Count; i++) { coeffs.a.Add(denominator[i] / scale); } for (int i = 0; i < numerator.Count; i++) { coeffs.b.Add(numerator[i] / scale); } return coeffs; }
All that needs to be done is scaling those values by gain, so let's add a function to Coefficients:
public class Coefficients { public List<double> a; public List<double> b; public void scaleGain(double gain) { for (int i = 0; i < b.Count; i++) { b[i] = b[i] / gain; } } }
Which gives us our completed lpf function in the new Butter filter.
// fcf is the cutoff frequency in Hz public void lpf(double fcf, int order) { double nfcf = fcf / PcmStreamSource.SampleRate; // normalised corner frequency Debug.Assert(nfcf > 0 && nfcf <= 0.5); // greater than 0 and less than Nyquist List<Complex> poles = getPoles(order); // get warped freq double w1 = 2 * Math.PI * Math.Tan(Math.PI * nfcf) / Math.PI; ComplexPlane splane = new ComplexPlane(); splane.poles = Polynomial.multiply(poles,w1); splane.zeros = new List<Complex>(); ComplexPlane zplane = splane.BilinearTransform(); TransferFunction transferfunction = zplane.getTransferFunction(); Coefficients coeffs = transferfunction.getCoefficients(); double gain = transferfunction.getGainAt(Complex.One); // DC Gain // normalise the gain by dividing coeffs.b by gain coeffs.scaleGain(gain); }
et voilà we have a new Butter filter which can specify low-pass functions.
So, let's take a look at the results. Generating some values for the following:
Butter b = new Butter(); b.lpf(10000, 2); b.lpf(10000, 3); b.lpf(10000, 4);
Plugging them into the Excel spreadsheet to graph the response, we get the following coefficients:
And a very nice looking graph:
No comments:
Post a Comment