Digital Filters in C#/Silverlight (Part 7 - ChebyshevII)

Hot on the heels of ChebyshevI, for completeness, let's also include the ChebyshevII filter. I'm not going to jump into maths again or even try to draw the geometric positioning of the poles and zeros today, so we'll just have to do with some code and graphs.

The simple bit is creating another class called Cheby2 and adding in the lpf function and getPoles using a suitable algorithm. Thankfully not too much work is required and Vincent Falco has already done a nice job in his Collection of Useful C++ Classes for DSP. In his code implementation he creates ComplexConjugate pairs and also Pole-Zero pairs, so we need to do a bit of additional work. We also need to modify the lpf function and add a new function in the Cheby2 class to generate some zeros as the ChebyshevII design requires them.

Just for completeness if browsing other code (including Vincent's), I've put the ChebyshevI algorithm in the more usual non-geometric calculation form below:

            // Cheby1,algorithm in more usual form
            for (int i = 0; i < order / 2; ++i)
            {
                int k = 2 * i + 1 - order;
                double x = sinhA * Math.Cos(k * Math.PI / (2.0 * (double)order));
                double y = coshA * Math.Sin(k * Math.PI / (2.0 * (double)order));

                Complex pole = new Complex(x, y);
                poles.Add(pole);
                poles.Add(Complex.Conjugate(pole));
            }

And here's the Cheby2 class:

using System;
using System.Diagnostics;
using System.Numerics;
using System.Collections.Generic;

// Hondrou implementation of a Chebyshev2 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
//
// The Chebyshev2 algorithm comes from the calculations by Vincent Falco in
// his Collection of Useful C++ Classes for DSP
// https://github.com/vinniefalco/DSPFilters 

namespace AgSynth
{
    public class Cheby2 : Iir
    {
        // ripple is stopband ripple in this case, but same calculation
        public List<Complex> getPoles(int order, double ripple)
        {
            Debug.Assert(ripple <= 0.0); // ripple must be negative

            // 10 is correct because we need the square
            double rip = Math.Pow(10.0, -ripple / 10.0); 
            double eps = Math.Sqrt(rip - 1.0);

            double a = asinh(1.0 / eps) / (double)order;

            Debug.Assert(a >= 0.0);

            double sinhA = -Math.Sinh(a);
            double coshA = Math.Cosh(a);

            List<Complex> poles = new List<Complex>();

            int k = 1;
            for (int i = order/2; --i >= 0; k+=2)
            {
                double x = sinhA * Math.Cos((k - order) * Math.PI / (2.0 * (double)order));
                double y = coshA * Math.Sin((k - order) * Math.PI / (2.0 * (double)order));

                double d2 = x * x + y * y;
                Complex pole = new Complex(x/d2, y/d2);
                poles.Add(pole);
                poles.Add(Complex.Conjugate(pole));
            }

            if (order % 2 != 0) //odd
            {
                poles.Add(Complex.FromPolarCoordinates(1 / sinhA, 0));
            }

            return poles;
        }

        public List<Complex> getZeros(int order, double ripple)
        {
            List<Complex> zeros = new List<Complex>();

            int k = 1;
            for (int i = order / 2; --i >= 0; k += 2)
            {
                double im = 1 / Math.Cos(k * Math.PI / (2.0 * (double)order));

                Complex zero = new Complex(0, im);
                zeros.Add(zero);
                zeros.Add(Complex.Conjugate(zero));
            }
            
            return zeros;
        }

        // fcf is the cutoff frequency in Hz
        public void lpf(double fcf, double ripple, 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, ripple);
            List<Complex> zeros = getZeros(order, ripple);

            // 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 = Polynomial.multiply(zeros, w1);

            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);

            setNumerator(coeffs.b);
            setDenominator(coeffs.a);
        }

    }
}

Generating some responses using this with the same -0.5 ripple as before for n=2,3,4

The ripple can be seen in the stopband.