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