Thursday 10 July 2014

Digital Filters in C#/Silverlight (Part 2 - some poles)

The pilot of a plane on its way out of Poland dies unexpectedly in flight. A passenger is asked to fill in. He looks at the controls and shakes his head. "What's wrong?" someone asks. The reply: "I'm just a simple Pole in a complex plane." 

Ahhh, the old ones are the best.....

Now we have the zeros under our belt, it's time to tackle their cousins, the poles. These are the ones that take the output and feedback scaled, delayed values. Let's start at the beginning with a single pole.

OnePole
The block diagram looks like this:








With the difference equation and transfer function:






With the following code:


using System;
using System.Net;

// Hondrou implementation of a One-Pole filter adapted from STK
// 
// This class implements a one-pole digital filter. A method is
// provided for setting the pole position along the real axis of the
// z-plane while maintaining a constant peak filter gain.
// by Perry R. Cook and Gary P. Scavone, 1995--2014.

namespace AgSynth
{
    public class OnePole : Filter
    {

        // The default constructor creates a low-pass filter (pole at z = 0.9).
        public OnePole()
        {
            b.Add(0.0); // b0
            a.Add(1.0); // a0
            a.Add(0.0); // a1

            setPole(0.9);

            inputs.Add(0.0); 
            outputs.Add(0.0); 
            outputs.Add(0.0);
        }

        // Set the pole position in the z-plane.
        // This method sets the pole position along the real-axis of the
        // z-plane and normalizes the coefficients for a maximum gain of one.
        // A positive pole value produces a low-pass filter, while a negative
        // pole value produces a high-pass filter. This method does not
        // affect the filter  gain value. The argument magnitude should be
        // less than one to maintain filter stability.
        public void setPole(double pole)
        {
            if (Math.Abs(pole) >= 1.0)
            {
                throw new Exception();
            }

            // normalise coefficients for peak unity gain
            if (pole > 0.0)
            {
                b[0] = (1.0 - pole);
            }
            else
            {
                b[0] = (1.0 + pole);
            }

            a[1] = -pole;
        }

        public double b0
        {
            set
            {
                b[0] = value;
            }
        }

        public double a1
        {
            set
            {
                a[1] = value;
            }
        }


        public void setCoefficients(double b0, double a1)
        {
            if (Math.Abs(a1) >= 1.0)
            {
                throw new Exception();
            }

            b[0] = b0;
            a[1] = a1;

        }

        public override double tick(double input)
        {

            inputs[0] = gain * input;
            last = b[0] * inputs[0] - a[1] * outputs[1];
            outputs[1] = last;

            return last;
        }
    }
}


It's worth looking at the tick function to see how the outputs are used and see that a[0] is never used in tick or the difference equation.

Moving swiftly on... Two poles...

A 747 was flying along and was full of Polish people. As they were going past some beautiful landmarks, the pilot came over the intercom and instructed all who were interested in seeing the landmark to look out the right side of the plane. Many passengers did so, and the plane promply crashed. Why? Too many poles in the right hand plane.

TwoPole
This has the following block diagram:











Difference equation and transfer function:

And code:


using System;
using System.Net;

// Hondrou implementation of a Two-Pole filter adapted from STK
// 
// This class implements a two-pole digital filter. A method is
// provided for creating a resonance in the frequency response while
// maintaining a nearly constant filter gain.
// by Perry R. Cook and Gary P. Scavone, 1995--2014.

namespace AgSynth
{
    public class TwoPole : Filter
    {
        // Default constructor creates a second-order pass-through filter.
        public TwoPole()
        {
            b.Add(1.0); // b0
            a.Add(1.0); // a0
            a.Add(0.0); // a1
            a.Add(0.0); // a2

            inputs.Add(0.0); 
            outputs.Add(0.0); 
            outputs.Add(0.0);
            outputs.Add(0.0);
        }


        public double b0
        {
            set
            {
                b[0] = value;
            }
        }

        public double a1
        {
            set
            {
                a[1] = value;
            }
        }

        public double a2
        {
            set
            {
                a[2] = value;
            }
        }

        public void setCoefficients(double b0, double a1, double a2)
        {
            if (Math.Abs(a1) >= 1.0)
            {
                throw new Exception();
            }

            b[0] = b0;
            a[1] = a1;
            a[2] = a2;

        }

        // Sets the filter coefficients for a resonance at frequency (in Hz).
        
        // This method determines the filter coefficients corresponding to
        // two complex-conjugate poles with the given frequency (in Hz)
        // and radius from the z-plane origin. If normalize is true,
        // the coefficients are then normalized to produce unity gain at
        // frequency (the actual maximum filter gain tends to be slightly
        // greater than unity when radius is not close to one). The
        // resulting filter frequency response has a resonance at the given
        // frequency. The closer the poles are to the unit-circle
        // (radius close to one), the narrower the resulting resonance width.
        // An unstable filter will result for \e radius >= 1.0. The
        // frequency value should be between zero and half the sample rate.
        // For a better resonance filter, use a BiQuad filter. 

        public void setResonance(double frequency, double radius)
        {
            if (frequency < 0.0 || frequency > 0.5 * PcmStreamSource.SampleRate)
            {
                throw new Exception();
            }

            if (radius < 0.0 || radius >= 1.0)
            {
                throw new Exception();
            }

            a[2] = radius * radius;
            a[1] = -2.0 * radius * Math.Cos(2 * Math.PI * frequency / PcmStreamSource.SampleRate);

            bool normalize = true;
            if (normalize)
            {
                // normalize the filter gain
                double real = 1.0 - radius + (a[2] - radius) * Math.Cos(Math.PI * 2.0 * frequency / PcmStreamSource.SampleRate);
                double imag = (a[2] - radius) * Math.Sin(Math.PI * 2.0 * frequency / PcmStreamSource.SampleRate);
                b[0] = Math.Sqrt(Math.Pow(real,2) + Math.Pow(imag,2));
            }
        }

        public override double tick(double input)
        {
            inputs[0] = gain * input;
            last = b[0] * inputs[0] - a[1] * outputs[1] - a[2] * outputs[2];
            outputs[2] = outputs[1];
            outputs[1] = last;

            return last;
        }
    }
}

The interesting property of adding poles is the possibility of generating some resonance through the feedback items. Using the function, again at 10kHz with some varying values, the following design paramters were created:
 
With the following Magnitude and Phase responses:


Which only leaves us with mixing both poles and zeros in the last foray into the STK digital filters.








No comments:

Post a Comment