Digital Filters in C#/Silverlight (Part 3 - poles and zeros)

Now the final bit of the story, mixing poles and zeros.

PoleZero
This has the following block diagram:


Difference equation and transfer function:

And code:

using System;
using System.Net;

// Hondrou implementation of a One-Pole filter adapted from STK
// 
// This class implements a one-pole, one-zero digital filter. A
// method is provided for creating an allpass filter with a given
// coefficient. Another method is provided to create a DC blocking
// filter.
// by Perry R. Cook and Gary P. Scavone, 1995--2014.

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

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


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

        public double b1
        {
            set
            {
                b[1] = value;
            }
        }

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

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

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

        // Set the filter for allpass behavior using coefficient.
        // This method uses coefficient to create an allpass filter,
        // which has unity gain at all frequencies. Note that the
        // coefficient magnitude must be less than one to maintain
        // filter stability.
        public void setAllpass(double coefficient)
        {
            if (Math.Abs(coefficient) >= 1.0)
            {
                throw new Exception();
            }

            b[0] = coefficient;
            b[1] = 1.0;
            a[0] = 1.0; // just in case
            a[1] = coefficient;
        }

          
        // Create a DC blocking filter with the given pole position in the z-plane.
        // This method sets the given pole position, together with a zero
        // at z=1, to create a DC blocking filter. The argument magnitude
        // should be close to (but less than) one to minimize low-frequency
        // attenuation.
        public void setBlockZero(double pole)
        {
            if (Math.Abs(pole) >= 1.0)
            {
                throw new Exception();
            }

            b[0] = 1.0;
            b[1] = -1.0;
            a[0] = 1.0; // just in case
            a[1] = -pole;
        }


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

            return last;
        }
    }
}

Now, let's have a look at some of the new interesting functions in polezero:

Allpass
Using the Excel sheet developed earlier we can plug in some of these values and take a look at the magnitude, phase and delay response:

This has a unitary gain magnitude response across all frequencies and the delay is linear(ish) over the frequency range. This is the transfer characteristic of the all-pass filter we previously saw in the DelayA class.

BlockZero
Again, using the Excel sheet we can take a look at the magnitude and phase response:

This attenuates, or blocks amplitude at and very close to zero frequency. This kind of filter is used to remove DC components from a signal.

Extending this the next obvious step is to add another pole and another zero making our very useful friend the BiQuad:

BiQuad
Which has the following block diagram:

Difference equation and transfer function:








And code:

using System;
using System.Net;

// Hondrou implementation of a Biquad filter adapted from STK
//
// This class implements a two-pole, two-zero digital filter.
// Methods are provided for creating a resonance or notch in the
// frequency response while maintaining a constant filter gain.
// by Perry R. Cook and Gary P. Scavone, 1995--2014.
namespace AgSynth
{
    public class BiQuad : Filter
    {
        // Default constructor creates a second-order pass-through filter.
        public BiQuad()
        {
            b.Add(1.0); // b0
            b.Add(0.0); // b1
            b.Add(0.0); // b2
            a.Add(1.0); // a0
            a.Add(0.0); // a1
            a.Add(0.0); // a2

            inputs.Add(0.0);
            inputs.Add(0.0);
            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 b1
        {
            set { b[1] = value; }
        }

        public double b2
        {
            set { b[2] = value; }
        }

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

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

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

            b[0] = b0;
            b[1] = b1;
            b[2] = b2;
            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 filter zeros are placed at z = 1, z = -1, and the coefficients
        // are then normalized to produce a constant unity peak gain
        // (independent of the filter gain parameter). 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 radius >= 1.0. The
        // frequency value should be between zero and half the sample rate.
        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.0 * Math.PI * frequency / PcmStreamSource.SampleRate);

            bool normalize = true;
            if (normalize)
            {
                // Use zeros at +- 1 and normalize the filter peak gain
                b[0] = 0.5 - 0.5 * a[2];
                b[1] = 0.0;
                b[2] = -b[0];
            }
        }

        // Set the filter coefficients for a notch at frequency (in Hz).
        // This method determines the filter coefficients corresponding to
        // two complex-conjugate zeros with the given frequency (in Hz)
        // and radius from the z-plane origin. No filter normalization is
        // attempted. The frequency value should be between zero and half
        // the sample rate. The radius value should be positive.
        public void setNotch(double frequency, double radius)
        {
            if (frequency < 0.0 || frequency > 0.5 * PcmStreamSource.SampleRate)
            {
                throw new Exception();
            }

            if (radius < 0.0)
            {
                throw new Exception();
            }

            // This method does not attempt to normalize the filter gain.
            b[2] = radius * radius;
            b[1] = -2.0 * radius * Math.Cos(2.0 * Math.PI * frequency / PcmStreamSource.SampleRate);
        }

        // Sets the filter zeroes for equal resonance gain.
        // When using the filter as a resonator, zeroes places at z = 1, z
        // = -1 will result in a constant gain at resonance of 1 / (1 - R),
        // where R is the pole radius setting.
        public void setEqualGainZeros()
        {
            b[0] = 1.0;
            b[1] = 0.0;
            b[2] = -1.0;
        }


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

            return last;
        }
    }
}

Biquad has the following interesting functions:

Resonance
Designs the filter coefficients to create a peak at the resonant frequency with the narrowness of the peak controlled by the radius paramters. It has the following response characteristics:

Finally bringing us to the general case, the Infinite Impulse Response (IIR) filter:

Iir
Which we've seen before, has the following block diagram:

Difference equation and transfer function:

 
And code:

using System;
using System.Net;
using System.Collections.Generic;

// Hondrou implementation of a Biquad filter adapted from STK
// This class provides a generic digital filter structure that can be
// used to implement IIR filters. For filters containing only
// feedforward terms, the Fir class is slightly more efficient.
//
// In particular, this class implements the standard difference
// equation:
//
// a[0]*y[n] = b[0]*x[n] + ... + b[nb]*x[n-nb] -
// a[1]*y[n-1] - ... - a[na]*y[n-na]
//
// If a[0] is not equal to 1, the filter coeffcients are normalized
// by a[0].
//
// The gain parameter is applied at the filter input and does not
// affect the coefficient values. The default gain value is 1.0.
// This structure results in one extra multiply per computed sample,
// but allows easy control of the overall filter gain.
//
// by Perry R. Cook and Gary P. Scavone, 1995--2014.

namespace AgSynth
{
    public class Iir : Filter
    {
        public Iir()
        {
            // default constructor sets up pass-through
            b.Add(1.0); // b0
            a.Add(1.0); // a0

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

        public Iir(List<double> bcoeff, List<double> acoeff)
        {
            b = bcoeff;
            a = acoeff;

            for (int i = 0; i < b.Count; i++)
            {
                inputs.Add(0.0);
            }

            for (int j = 0; j < a.Count; j++)
            {
                outputs.Add(0.0);
            }

            gain = 1.0;
        }


        public void setCoefficients(List<double> bcoeff, List<double> acoeff)
        {
            setNumerator(bcoeff);
            setDenominator(acoeff);
        }

        public void setNumerator(List<double> bcoeff)
        {
            b = bcoeff;

            inputs.Clear();
            for (int i = 0; i < b.Count; i++)
            {
                inputs.Add(0.0);
            }
        }

        public void setDenominator(List<double> acoeff)
        {
            a = acoeff;

            outputs.Clear();
            for (int i = 0; i < a.Count; i++)
            {
                outputs.Add(0.0);
            }

            // scale coefficients by a[0] if necessary
            if (a[0] != 1.0)
            {
                for (int i = 0; i < b.Count; i++)
                {
                    b[i] /= a[0];
                }
                for (int i = 1; i < a.Count; i++)
                {
                    a[i] /= a[0];
                }
            }
        }


        public override double tick(double input)
        {
            outputs[0] = 0.0;
            inputs[0] = input * gain;

            for (int i = b.Count - 1; i > 0; i--)
            {
                outputs[0] += b[i] * inputs[i];
                inputs[i] = inputs[i - 1];
            }

            outputs[0] += b[0] * inputs[0];

            for (int i = a.Count - 1; i > 0; i--)
            {
                outputs[0] += -a[i] * outputs[i];
                outputs[i] = outputs[i - 1];
            }

            last = outputs[0];

            return last;
        }
    }
}