Thursday 24 July 2014

Digital Filters in C#/Silverlight (Part 6 - ChebyshevI)

Having got the Butterworth filter under our belt it makes obvious sense to get on cracking and implement a ChebyshevI filter design as well as the design technique follows on pretty much from the established groundwork. We're not going to lurch back into theory except to explain how the ChebyshevI design comes about.
As we saw earlier, the Butterworth design places poles onto the s-plane around a unit circle in an ordered pattern. The ChebyshevI design places these on an ellipse in a similar configuration which produces a steeper cutoff characteristic but introduces ripples into the passband. The amount of ripple is adjusted by the 'squashing' of the ellipse towards the imaginary axis and it's obvious to see that a Butterworth design therefore is simply a special case of the ChebyshevI design where the ellipse is a perfect circle (e.g. the ripple parameter = 0).

The s-plane placement therefore looks something like this:













Where \( \alpha_{max} \) is the ripple and:

$$ a = \frac {sinh^{-1}\left(\frac {1}{\epsilon}\right)}{n} $$
$$ \epsilon = \sqrt {10^\frac {\alpha_{max}}{10} - 1} $$

The code makes use of  Guillemin's algorithm for pole placement which shifts the pole positions geometrically from the Butterworth positions rather than calculating them directly from the Chebyshev polynomial.

So, this makes the coding pretty easy. The lpf, bpf, bpf, bsf functions can all remain pretty much the same except we need to create a new getPoles function that has a parameter for ripple. I'm not going to show the whole class as it should be pretty much cut-and-paste, but here's the new getPoles function for Cheby1:


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

// Hondrou implementation of a Chebyshev1 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 Cheby1 : Iir
    {
        private static double asinh(double x)
        {
            return Math.Log(x + Math.Sqrt(1.0 + x * x));
        }

        public List<Complex> getPoles(int order, double ripple)
        {
            Debug.Assert(ripple <= 0.0); // ripple must be negative

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

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

            Debug.Assert(y >= 0.0);

            double sinhY = Math.Sinh(a);
            double coshY = Math.Cosh(a);

            Complex[] poles = new Complex[order];

            for (int i = 0; i < order; i++)
            {
                double theta = (order / 2.0 + 0.5 + i) * Math.PI / order;
                Complex circle = Complex.FromPolarCoordinates(1, theta);
                poles[i] = new Complex(circle.Real * sinhY, circle.Imaginary * coshY);
            }

            return new List<Complex>(poles);
        }

        // the rest of the lpf, hpf, bpf, bsf functions go here
    }
}

And here are the example responses using the same values previously with the Butterworth design and ripple set to -0.5.

Lowpass (lpf)
















Bandpass (bpf)
















Bandstop (bsf)




No comments:

Post a Comment