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

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.

Digital Filters in C#/Silverlight (Part 1 - a bunch of Zeros)

As I mentioned previously last night I had a blitz and implemented all of the filters in STK so that I now have a range of audio tools to play with and wanted to get this down in the blog. I'm going to break this into a number of chunks to walk through the implementation and some of the concepts.

Kicking off, we'll start from the basics and build in more functionality, so the best place to start with is the OneZero filter, named because it has one zero!

OneZero
This filter is our old friend the two point averager, it takes the current input, scales it and adds it to a scaled delayed input. Here's the block diagram:

And the difference equation and transfer function:

The code looks like this:

using System;
using System.Net;

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

namespace AgSynth
{
    public class OneZero : Filter
    {
        // default zero at -1.0
        public OneZero()
        {
            b.Add(0.0); // b0
            b.Add(0.0); // b1

            setZero(-1.0);

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

        // Set the zero position in the z-plane.
        // This method sets the zero position along the real-axis of the
        // z-plane and normalizes the coefficients for a maximum gain of one.
        // A positive zero value produces a high-pass filter, while a
        // negative zero value produces a low-pass filter. This method does
        // not affect the filter gain value.
        public void setZero(double zero)
        {
            // normalise coefficients for peak unity gain
            if (zero > 0.0)
            {
                b[0] = 1.0 / (1.0 + zero);
            }
            else
            {
                b[0] = 1.0 / (1.0 - zero);
            }

            b[1] = -zero * b[0];
        }

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

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


        public void setCoefficients(double b0, double b1)
        {
            b[0] = b0;
            b[1] = b1;
        }

        public override double tick(double input)
        {

            inputs[0] = gain * input;
            last = b[1] * inputs[1] + b[0] * inputs[0];
            inputs[1] = inputs[0];

            return last;
        }
    }
}

See how the tick operation mirrors the difference equation.

Moving swiftly on, you guessed it, yep, TwoZeros:

TwoZero
You can guess the pattern, by now. Here is the block diagram:

And the difference equation and transfer function:


And code:

using System;
using System.Net;

// Hondrou implementation of a Two-Zero filter adapted from STK
//
// This class implements a two-zero digital filter. A method is
// provided for creating a "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 TwoZero : Filter
    {
        // Default constructor creates a second-order pass-through filter.
        public TwoZero()
        {
            b.Add(1.0); // b0
            b.Add(0.0); // b1
            b.Add(0.0); // b2

            inputs.Add(0.0); 
            inputs.Add(0.0);
            inputs.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 void setCoefficients(double b0, double b1, double b2)
        {
            b[0] = b0;
            b[1] = b1;
            b[2] = b2;
        }

        // Sets 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. The coefficients are then
        // normalized to produce a maximum filter gain of one (independent of
        // the filter gain parameter). The resulting filter frequency
        // response has a "notch" or anti-resonance at the given
        // frequency. The closer the zeros are to the unit-circle (radius
        // close to or equal to one), the narrower the resulting notch width.
        // 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();
            }

            b[2] = radius * radius;
            b[1] = -2.0 * radius * Math.Cos(2.0 * Math.PI * frequency / PcmStreamSource.SampleRate);
        
            // Normalise the filter gain
            if (b[1] > 0.0) // max at z=0
            {
                b[0] = 1.0 / (1.0 + b[1] + b[2]);
            }
            else
            {
                b[0] = 1.0 / (1.0 - b[1] + b[2]);
            }
            b[1] *= b[0];
            b[2] *= b[0];
        }


        public override double tick(double input)
        {

            inputs[0] = gain * input;
            last = b[2] * inputs[2] + b[1] * inputs[1] + b[0] * inputs[0];
            inputs[2] = inputs[1];
            inputs[1] = inputs[0];

            return last;
        }
    }
}

What is interesting about this filter implementation is the notch function.

Notch
Using the Excel sheet from my previous post with some tweaks to show multiple filters, we can now look at how the filter can be designed to perform in certain ways. I've used a notch frequency of 10,000 Hz and set the radius value to some different values.

The filter design parameters look like this:

 

Giving the following Magnitude and Phase response:

Pretty useful!

This just leaves us with the general case, our friend the FIR:

Fir
For completeness, the block diagram:


Difference equation and transfer function



And code

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

// Hondrou implementation of a Fir filter adapted from STK
//
// This class provides a generic digital filter structure that can be
// used to implement FIR filters. For filters with feedback terms,
// the Iir class should be used.
// In particular, this class implements the standard difference
// equation:
//
// y[n] = b[0]*x[n] + ... + b[nb]*x[n-nb]
//
// 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 Fir : Filter
    {
        // The default constructor should setup for pass-through.
        public Fir()
        {
            // default constructor sets up pass-through
            b.Add(1.0); // b0

            inputs.Add(0.0); 
        }

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

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

            gain = 1.0;
        }


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

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


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

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

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

            return last;
        }
    }
}

Once more unto the Averager (and IIR responses in Excel)

Having teased that I'd get round to explaining the Averager frequency response in my last blog, we're actually going to get to it now. There was a little more groundwork to cover than I expected and it was better to put into two posts.

We now have the tools in the Filter class to be able to do the phaseDelay and magnitude response for a filter. Just to make things a little easier, rather than putting a lot of additional code into our libraries I've taken the same algorithms and dropped them into Excel. It's a nice way of sharing the calculation algorithms again and actually is a neat little reference tool that can be used without running up a compiler.

When I have a spare evening I'll make a little app and put the graphs online in Silverlight so that this analysis can be used in some later things we're going to look at. For now, Excel it is, which can be downloaded here.

The spreadsheet allows the a and b variables for an IIR to be set in the green coloured columns.

Here you can see that b0 and b1 have been set to 0.5, which is our two-point averager and the number of coefficients in the filter can be set in the pink boxes, in this case just two delay terms. The number of feedback terms has been set to zero.











Finally, to drum-roll, the magnitude, phase and delay curves look like this:

So, you can see that the magnitude (amplitude) response does indeed show a low-pass characteristic, although not with much of a roll-off. The phase response is linear and the phase delay is constant at 0.5 (as we had stated).

A note on the frequency values. The scale is shown from 0-22.1kHz which is the Nyquist limit for a 44100 samples per second digital system used in audio. The phase response shows in this case just negative, with the scale from +pi to -pi in radians.

Now we have the code and a spreadsheet to go forth and start playing with filters.

Back to the Averager (and a bit of digital filter groundwork for lunch)

Well, last night I wanted to noodle around with some more sounds and get the basic KS on a little further so ended up deep in digital filter land, more of which later. I ended up the evening having implemented all the digital filter classes in the STK as they proved to be a good set of building blocks and an easy way to bring the concepts together in an hour or so. I'm going to blog on that next, but before going there I wanted to tie-up the lose end I mentioned with the Averager and in doing so take the opportunity to dip into digital filter land. So, this is kind of an explanation and pre-step to the next journey.

Intuitively it's obvious that a basic two point moving averager is a low-pass filter, and we did the hand-waving bit by putting the phase delay to be 0.5 without taking that any further. Before going into digital filters proper, it's worth understanding these two points (ha ha) as the thinking is good to explore the more complicated stuff later on.

As an online introduction reference for digital filters, keeping with the STK and audio theme, the best reference is Julius O Smith's Introduction to Digital Filters with Audio Applications.

So, stepping straight in, we're going to go for the general case here to setup some tools and then when I blog about the STK filter implementation we'll use those tools and build upwards. Constraining the problem, we're only going to look at Linear Time Invariant (LTI) filters. The general form we'll be looking at are Infinite Impulse Response (IIR) filters.

IIR filter
In block diagram terms they look like the diagram below. In simple terms there are a set of delayed inputs, each scaled by a value b added together with a set of delayed outputs (feedback) each scaled by a value a.

What this means is many of the different forms of delay and filter (a delay is just a kind of filter) we have already seen and will go on to look at can all be generalised to the same model. A simplification of this model where there is no feedback present and all a values are set to zero is called a Finite Impulse Response (FIR) filter. For example the simple non-interplotating delay line could be considered as a FIR filter where the delay items b0,b1,.... b(N-1) are all zero and bN is given the value 1. This would just delay the input by the number of samples N.

Back to our two point moving averager, this is actually a FIR filter with b0 = 0.5 and b1 = 0.5. Explaining this, b0 is the current input and we take the last input (delayed by one sample) add the two together and divide by two (averaging) or as we can see, mix half of each together.

So, in our Averager class implementation as this is a subclass of Filter we can set the two parameters b[0] and b[1] = 0.5 which will allow us to use some of the other tools we're going to look at. It will also mean that the overriden phaseDelay function can be removed and the general one in Filter can be used instead as we'll see a bit later on.

        public Averager()
        {
            b = new List<double>() { 0.5, 0.5 };
        }

Before going further, we should have a look at how the block diagram seems in notation which will help us navigate some of the maths going on and the code algorithms after that. The difference equation and transfer function for an IIR are shown below:

As there are a lot of symbols here to start with, I've taken the small simplification of making a0=1 which is the usual case as it looks slightly neater and will not get in the way where we're going. Please consult one of the references to dig deeper if needed.

What the equations are saying is that the output is a sum of the scaled (by b-values) of delayed inputs x(n-p) along with scaled (by a-values) of the feedback outputs y[n-q].

So if you set all the a's to zero half of the equation disappears like the block diagram. 

Frequency Response

Now, what we're interested in with filters is getting some kind of frequency response. What this means is that for certain frequencies we would like the input signal to have the amplitude of those frequencies scaled by a certain amount. I'm assuming here that this grounding of that is in place, so will not distract us by going into details in this post. Now, this is the usual departure point to go into lots of maths, but I'll leave that for the textbooks. What we want is some way numerically (e.g. in code) of being able to plug a set of these a and b parameters into a model and look at how the filter affects certain frequencies to see what effect we're getting.

Robin Schmit has a nice little document on his site with some of the equations we need. The one below calculates the magnitude response of a filter at a certain frequency:

This is basically taking apart bits of the equations and considering the phasors at a particular frequency and seeing how they are scaled. I'm not going to go any deeper with this as the equation is enough to get us back into code.

The other companion equation is the phase response.

This is like the magnitude response, but shows how the phase of a certain input frequency signal is affected by the filter.

Now, bringining this together. Remember that equation phaseDelay that we glossed over in the Filter class? Hmmm, well, looks pretty similar to the phase response, right? Nope, well take a look, it's summing up the reala and imaga (which is ca and sa in the above equations) and then doing a little atan2 calculation.

Well, we're nearly there, the eagle eyed will have noticed that this equation is the phase response and the code is calculating the phaseDelay. Well, that's simply just a shift in terms to see how much delay the phase shift means, the concepts are interchangable. 

So, we'd better get back into some code and put in this magnitude response function, which can be added to Filter:

        public virtual double magnitude(double frequency)
        {
            if (frequency <= 0.0 || frequency > 0.5 * PcmStreamSource.SampleRate)
            {
                throw new Exception();
            }

            double omegaT = 2 * Math.PI * frequency / PcmStreamSource.SampleRate;
            double reala = 1.0;
            double imaga = 0.0;
            double realb = 0.0;
            double imagb = 0.0;
 
            for (int i = 0; i < a.Count; i++)
            {
                reala += a[i] * Math.Cos(i * omegaT);
                imaga -= a[i] * Math.Sin(i * omegaT);
            }

            for (int i = 0; i < b.Count; i++)
            {
                realb += b[i] * Math.Cos(i * omegaT);
                imagb -= b[i] * Math.Sin(i * omegaT);
            }

            reala *= gain;
            realb *= gain;
            imaga *= gain;
            imagb *= gain;

            double mag = Math.Sqrt((Math.Pow(realb,2) + Math.Pow(imagb,2)) / 
                                   (Math.Pow(reala,2) + Math.Pow(imaga,2)));

            return mag;
        }

With this in place we're ready to rock and roll and take a look at the magnitude and phase response for the averager and any of these type of filters that can be made from this in the future.

It's a bit kludgy, but here's a simple algorithm to show using the functions with the averager and putting the values into a long text string. Why? Well it's an easy way to get the values out and into some graphing tools.

         private void AnalyseResponse()
        {
            Averager filter = new Averager();

            string phase = "";
            string amp = "";

            int SAMPLES = 100;
            for (int i = 1; i < SAMPLES; i++)
            {
                double ph = filter.phaseDelay(PcmStreamSource.SampleRate / 2.0 * ((double)i / (double)SAMPLES));
                double a = filter.amplitude(PcmStreamSource.SampleRate / 2.0 * ((double)i / (double)SAMPLES));

                phase += String.Format("{0:0.00},", ph);
                amp += String.Format("{0:0.00},", a);
            }
        }

Having promised to explain the averager response, I'm not quite going to get there in this post. So, finishing up for lunch, I'll quickly post how to show this response in another blog note.

Bringing the Generator classes closer to STK

After the work just done on the Delay classes the effort can now really build some benefit by harmonising the Generator classes more closely with STK and the tick() method. We'll do this progressively and then look at making some further changes in a follow-on post.

First of all, until now the tick method has been returning a sample value to the PcmStreamSource as a short value. We're going to change it to return a double and then put the scaling into PcmStreamSource. This keeps all the generator code in the floating point domain.

    public interface ISampleSource
    {
        double tick();
    }

and then in GetSampleAsync (not all the function is shown for brevity:

        protected override void GetSampleAsync(MediaStreamType mediaStreamType)
        {
            for (int i = 0; i < numSamples; i++)
            {
                short sample = 0;
                if (this.Input != null)
                {
                    sample = (short)(this.Input.tick() * short.MaxValue);
                }

                for (int c = 0; c < ChannelCount; c++)
                {
                    //channel
                    memoryStream.WriteByte((byte)(sample & 0xFF));
                    memoryStream.WriteByte((byte)(sample >> 8));
                }
            }

The corresponding changes then need to be made in the Generator interface and in the Karp implementation:

    public class Generator : ISampleSource
    {
        public virtual double tick()
        {
            return 0.0;
        }
    }

The changes to the Karp class are quite obvious to change this to a double so I will not explicitly describe them as we've got some more fundamental changes on the way. What we're going to do is unpack the various different filtering components now from the Karp class and put them into separate filters and generators so we have more room to create a wider range of generators in the future like the STK and develop a wider range of filters.

Noise.cs
First of all, let's look at making a new Noise generator class and take the functionality outside of the Karp class:

using System;
using System.Net;


namespace AgSynth
{
    public class Noise : Generator
    {
        private int noise = 0; // gate duration for generating noise

        Random r = new Random();


        // duration to generate noise in samples
        public int Duration
        {
            set
            {
                noise = value;
            }

            get
            {
                return noise;
            }
        }

        public override double tick()
        {
            double y = 0;

            // push noise into circuit if still generating noise
            if (noise > 0)
            {
                y = (r.NextDouble() * 2.0) - 1.0; // generate noise sample
                noise--;
            }

            return y;
        }
    }
}

This is pretty simple, but slightly different from the STK. In the STK Plucked implementation, which is closest to our basic Karp class noise is added into the delay when plucked is called. We might come back and take a look at the pros and cons of this, for now we'll stick with our current implemtation which effectively generates a duration of noise which is provided every tick.

Averager.cs
One of the reasons for unpacking the functions to match more closely to STK is so that we can have a look at a wider range of filters to build up some more interesting instruments. For now, rather than complicate things let's make a new filter just for the two-point moving averager which is used for the LPF. This is subclass of Filter. For simplicity as well we're setting the phaseDelay explicitly to 0.5 rather than using the calculation. I'll try to cover that in a separate post later on. It's more important that we take a look at how the Karp class is going to change. For now the Averager looks like this:

using System;
using System.Net;


namespace AgSynth
{
    public class Averager : Filter
    {
        public override double phaseDelay(double frequency)
        {
            return 0.5;
        }

        private double ydel = 0.0;

        public override double tick(double input)
        {
            double y = (input + ydel) * 0.5 * gain;
            ydel = y;

            return y;
        }
    }
}

I'm trying to keep the notation in the tick functions close to using the x and y values from the difference equations rather than other variables so that the theory is easier to follow.

Now the changes to Karp, which becomes much leaner. A lot of the variables are teased out and the two new classes Noise and Averager are added as attributes. The gain of the Averager (lpf) is set in the constructor to the 0.9999 value used previously. Pluck is now used to kick-off the call to set the noise generation duration using the delay.Delay length as a variable. The hard-coded value for the lpf phase delay in the calculation of the length of the delay line has now using the Filter phaseDelay function which will help us in the future. The biggest changes are to the tick operation which now unbundles the algorithm into the various generators and filters. I'm going to show the whole class so all the changes can be seen and then we'll look at a final little change.

using System;
using System.Net;

namespace AgSynth
{
    public class Karp: Generator
    {
        double Fo = 0; // fundamental frequency (pitch requested)

        DelayA delay = new DelayA();
        Noise noise = new Noise();
        Averager lpf = new Averager();

        public Karp()
        {
            lpf.Gain = 0.9999;
        }

        public void pluck()
        {
            noise.Duration += Helpers.Truncate(delay.Delay * 0.95); ;
        }


        public void pitch(double p)
        {
            if (p == 0)
            {
                delay.Delay = 0;
            }
            else
            {
                Fo = p; // frequency
                double Fs = PcmStreamSource.SampleRate;

                double No = Fs / Fo;

                delay.Delay = No - lpf.phaseDelay(p);
            }
        }



        public override double tick()
        {
            double sample = 0;

            if (delay.Delay > 0)
            {
                double amp = 0.99; // amplification

                // get input which is noise added with feedback loop
                double x = delay.pop() + noise.tick();  

                // filter the output with simple LPF
                double y = lpf.tick(x);

                // push the output to the start of the delay-line
                delay.tick(y);

                // amplify
                sample = (y * amp);
            }

            return sample;
        }

    }
}

In the STL the tick operation is a little more bundled. I've kept with the x and y notation as I think this is a little easier to read and understand.

The final little bit is that we can move a Gain value into the Generatorclass which can be used in Noise and Karp.

    public class Generator : ISampleSource
    {
        protected double gain = 1.0;

        public double Gain
        {
            set
            {
                gain = value;
            }
            get
            {
                return gain;
            }
        }

        public virtual double tick()
        {
            return 0.0;
        }
    }

The tick operation in Noise can now be changed to add in gain:

        public override double tick()
        {
            double y = 0;

            // push noise into circuit if still generating noise
            if (noise > 0)
            {
                y = (r.NextDouble() * 2.0 * gain) - 1.0; // generate noise sample
                noise--;
            }

            return y;
        }

This allows the pluck operation in Karp to be given an amplitude value (how hard the string is plucked) by adjusting the gain of the noise generator. We can see how this can also be used for some different effects later.

        public void pluck(double a)
        {
            noise.Gain = a;
            noise.Duration += Helpers.Truncate(delay.Delay * 0.95); ;
        }

And likewise the gain value can be used for the Karp tick operation itself:

        public override double tick()
        {
            double sample = 0;

            if (delay.Delay > 0)
            {
                // get input which is noise added with feedback loop
                double x = delay.pop() + noise.tick();  

                // filter the output with simple LPF
                double y = lpf.tick(x);

                // push the output to the start of the delay-line
                delay.tick(y);

                // amplify
                sample = y * gain;
            }

            return sample;
        }

So, this refactoring gives no real difference in what is being generated (so no new app to show), but it builds a stronger framework that we can now explore and reference the STK a bit more closely. More in due course....