Digital Filters in C#/Silverlight (Part 4 - RBJ Cookbook)

Having completed implementing the STK filters, to use some of these in anger we now need to be able to set the coefficients to put them to any meaningful use. To the ready, the oft cited Cookbook formulae for audio EQbiquad filter coefficients written by Robert Bristow-Johnson.

I've implemented a sub-class from BiQuad called Rbj that has a set of functions to set the coefficients for the various filter forms. I'm going to go through the different functions first and then put the code in at the end.

Low-Pass / High-Pass Filter

Using the formulae from the cookbook and puting them into the Excel spreadsheet, the following kind of functions can be created.

 
This is showing low-pass over a different set of frequencies. In each case the Q value is set to 0.9. We can now take a look at what changing the Q value can do for a single frequence (10k).

In this case, setting the Q value to 1/SQRT(2) makes the function the standard Butterworth rolloff. Values above this show a resonant peak which can be interesting for some applications.

 

The high-pass function is effectively a symmetric mirror of this.

Bandpass / Notch

The usual bandpass and notch configurations are also provided and look like this:

With varying frequency cut-offs changing as follows:

And the corresponding notch functions looking like this:






With the following frequency changes:

Coefficients

The equivalent analog transfer function and the coefficient paramters are shown below:

with the following paramters:

where f0 is the corner frequency of interest and fs is the sampling frequency 44100Hz in our case.

where Q is the quality factor for adjusting the resonance (or inverse damping) or in the case of the bandpass/notch can be defined in terms of the bandwidth BW.

don't forget to normalise a0 as needed!

Code

And finally the code looks like this:

using System;
using System.Net;
using System.Diagnostics;

// Hondrou implementation of a Rbj filters adapted from cookbook
//
// Cookbook filters for audio EQ biquad filter coefficients
//
// Implementation of BiQuad functions with coefficients from 
// Robert Bristow-Johnson  <rbj@audioimagination.com>
//
// All filter transfer functions were derived from analog prototypes (that
// are shown below for each EQ filter type) and had been digitized using the
// Bilinear Transform.  BLT frequency warping has been taken into account for
// both significant frequency relocation (this is the normal "prewarping" that
// is necessary when using the BLT) and for bandwidth readjustment (since the
// bandwidth is compressed when mapped from analog to digital using the BLT).
//
// Q (the EE kind of definition, except for peakingEQ in which A*Q is
// the classic EE Q.  That adjustment in definition was made so that
// a boost of N dB followed by a cut of N dB for identical Q and
// f0/Fs results in a precisely flat unity gain filter or "wire".)

namespace AgSynth
{
    public class Rbj : BiQuad
    {
        private void Normalise()
        {
            b[0] /= a[0];
            b[1] /= a[0];
            b[2] /= a[0];
            a[1] /= a[0];
            a[2] /= a[0];
        }

        // H(s) = 1 / (s^2 + s/Q + 1)
        public void setLpf(double f0, double Q)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = (1.0 - Math.Cos(w0)) / 2.0;
            b[1] = 1.0 - Math.Cos(w0);
            b[2] = (1.0 - Math.Cos(w0)) / 2.0;
            a[0] = 1.0 + alpha;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 1.0 - alpha;

            Normalise();
        }

        // H(s) = s^2 / (s^2 + s/Q + 1)
        public void setHpf(double f0, double Q)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = (1.0 + Math.Cos(w0)) / 2.0;
            b[1] = -(1.0 + Math.Cos(w0));
            b[2] = (1.0 + Math.Cos(w0)) / 2.0;
            a[0] = 1.0 + alpha;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 0.0;

            Normalise();
        }

        // H(s) = s / (s^2 + s/Q + 1)  (constant skirt gain, peak gain = Q)
        public void setBpf1(double f0, double Q)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = Math.Sin(w0) / 2.0; // = Q * alpha
            b[1] = 0.0;
            b[2] = -Math.Sin(w0) / 2.0; // -Q * alpha
            a[0] = 1.0 + alpha;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 1 - alpha;

            Normalise();
        }


        // H(s) = (s/Q) / (s^2 + s/Q + 1)      (constant 0 dB peak gain)
        public void setBpf2(double f0, double Q)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = alpha;
            b[1] = 0.0;
            b[2] = -alpha;
            a[0] = 1.0 + alpha;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 1 - alpha;

            Normalise();
        }

        // H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
        public void setNotch(double f0, double Q)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = 1.0;
            b[1] = -2.0 * Math.Cos(w0);
            b[2] = 1.0;
            a[0] = 1.0 + alpha;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 1 - alpha;

            Normalise();
        }

        // H(s) = (s^2 - s/Q + 1) / (s^2 + s/Q + 1)
        public void setApf(double f0, double Q)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = 1.0 - alpha;
            b[1] = -2.0 * Math.Cos(w0);
            b[2] = 1.0 + alpha;
            a[0] = 1.0 + alpha;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 1 - alpha;

            Normalise();
        }


        // H(s) = (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1)
        public void setPeakingEQ(double f0, double Q, double A)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = 1.0 + alpha * A;
            b[1] = -2.0 * Math.Cos(w0);
            b[2] = 1.0 - alpha * A;
            a[0] = 1.0 + alpha / A;
            a[1] = -2.0 * Math.Cos(w0);
            a[2] = 1 - alpha / A;

            Normalise();
        }


        // H(s) = A * (s^2 + (sqrt(A)/Q)*s + A)/(A*s^2 + (sqrt(A)/Q)*s + 1)
        public void setLowShelf(double f0, double Q, double A)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = A * ((A+1.0) - (A-1.0)*Math.Cos(w0) + 2.0 * Math.Sqrt(A)*alpha);
            b[1] = 2.0*A*((A+1.0) - (A+1.0)*Math.Cos(w0));
            b[2] = A * ((A+1.0) - (A-1.0)*Math.Cos(w0) - 2.0 * Math.Sqrt(A)*alpha);
            a[0] = (A+1.0) + (A-1.0)*Math.Cos(w0)+2.0*Math.Sqrt(A)*alpha;
            a[1] = -2.0 * ((A+1.0)+(A+1.0)*Math.Cos(w0));
            a[2] = (A+1.0) + (A-1.0)*Math.Cos(w0) - 2.0*Math.Sqrt(A)*alpha;

            Normalise();
        }

        // H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)
        public void setHighShelf(double f0, double Q, double A)
        {
            double w0 = 2.0 * Math.PI * f0 / PcmStreamSource.SampleRate;
            double alpha = Math.Sin(w0) / (2.0 * Q);

            b[0] = A * ((A+1.0) + (A-1.0)*Math.Cos(w0) + 2.0 * Math.Sqrt(A)*alpha);
            b[1] = -2.0*A*((A+1.0) + (A+1.0)*Math.Cos(w0));
            b[2] = A * ((A+1.0) + (A-1.0)*Math.Cos(w0) - 2.0 * Math.Sqrt(A)*alpha);
            a[0] = (A+1.0) - (A-1.0)*Math.Cos(w0)+2.0*Math.Sqrt(A)*alpha;
            a[1] = -2.0 * ((A+1.0) - (A+1.0)*Math.Cos(w0));
            a[2] = (A+1.0) - (A-1.0)*Math.Cos(w0) - 2.0*Math.Sqrt(A)*alpha;

            Normalise();
        }



    }
}