Lowpass, Highpass, and Bandpass Butterworth Filters in C#

• Download Butterworth_Filter_Design.zip - 17.5 KB

Introduction

This is a C# implementation of digital lowpass, highpass, and bandpass Butterworth filters of arbitrary order (n cascaded 2-pole sections).

Background

A Word document giving the filter design via bilinear z-transformation is included.

The design starts with a continuous-time lowpass Butterworth filter and uses the bilinear z-transform to derive the coefficients of the equivalent digital filter. Then a substitution of variable s=1/s transforms the lowpass filter into a highpass filter, and the bilinear z-transform is used to derive the coefficients of the equivalent highpass digital filter.

In the code, the lowpass and highpass filters are implemented according to the coefficient derivation from the design document. A bandpass filter is implemented as a cascade of a lowpass and a highpass filter.

The code includes an implementation of an n^th-order FIR filter for the zero (numerator) polynomials and an implementation of an n^th-order IIR filter for the pole (denominator) polynomials.

Using the Code

The code is divided into five files:

1. LowpassFilterButterworthImplementation.cs
2. HighpassFilterButterworthImplementation.cs
3. BandpassFilterButterworthImplementation.cs
4. FIRFilterImplementation.cs
5. IIRFilterImplementation.cs

To use the LowpassFilterButterworthImplementation, HighpassFilterButterworthImplementation, or BandpassFilterButterworthImplementation, construct the object with the desired cutoff frequency or frequencies, number of 2-pole sections, and the sample frequency (Fs). Then call the compute() method on the object repeatedly, passing in an input sample and receiving a filtered sample as the return value.

using DSP;
using System;

namespace DSP
{
    public class LowpassFilterButterworthImplementation
    {
        protected LowpassFilterButterworthSection[] section;

        public LowpassFilterButterworthImplementation
               (double cutoffFrequencyHz, int numSections, double Fs)
        {
            this.section = new LowpassFilterButterworthSection[numSections];
            for (int i = 0; i < numSections; i++)
            {
                this.section[i] = new LowpassFilterButterworthSection
                                  (cutoffFrequencyHz, i + 1, numSections * 2, Fs);
            }
        }
        public double compute(double input)
        {
            double output = input;
            for (int i = 0; i < this.section.Length; i++)
            {
                output = this.section[i].compute(output);
            }
            return output;
        }
    }
}

public class LowpassFilterButterworthSection
{
    protected FIRFilterImplementation firFilter = new FIRFilterImplementation(3);
    protected IIRFilterImplementation iirFilter = new IIRFilterImplementation(2);

    protected double[] a = new double[3];
    protected double[] b = new double[2];
    protected double gain;

    public LowpassFilterButterworthSection
           (double cutoffFrequencyHz, double k, double n, double Fs)
    {
        // compute the fixed filter coefficients
        double omegac = 2.0 * Fs * Math.Tan(Math.PI * cutoffFrequencyHz / Fs);
        double zeta = -Math.Cos(Math.PI * (2.0 * k + n - 1.0) / (2.0 * n));

        // fir section
        this.a[0] = omegac * omegac;
        this.a[1] = 2.0 * omegac * omegac;
        this.a[2] = omegac * omegac;

        //iir section
        //normalize coefficients so that b0 = 1, 
        //and higher-order coefficients are scaled and negated
        double b0 = (4.0 * Fs * Fs) + (4.0 * Fs * zeta * omegac) + (omegac * omegac);
        this.b[0] = ((2.0 * omegac * omegac) - (8.0 * Fs * Fs)) / (-b0);
        this.b[1] = ((4.0 * Fs * Fs) - 
                     (4.0 * Fs * zeta * omegac) + (omegac * omegac)) / (-b0);
        this.gain = 1.0 / b0;
    }

    public double compute(double input)
    {
        // compute the result as the cascade of the fir and iir filters
        return this.iirFilter.compute
               (this.firFilter.compute(this.gain * input, this.a), this.b);
    }
}

using DSP;
using System;

namespace DSP
{
    public class HighpassFilterButterworthImplementation
    {
        protected HighpassFilterButterworthSection[] section;

        public HighpassFilterButterworthImplementation
        (double cutoffFrequencyHz, int numSections, double Fs)
        {
            this.section = new HighpassFilterButterworthSection[numSections];
            for (int i = 0; i < numSections; i++)
            {
                this.section[i] = new HighpassFilterButterworthSection
                (cutoffFrequencyHz, i + 1, numSections * 2, Fs);
            }
        }
        public double compute(double input)
        {
            double output = input;
            for (int i = 0; i < this.section.Length; i++)
            {
                output = this.section[i].compute(output);
            }
            return output;
        }
    }
}

public class HighpassFilterButterworthSection
{
    protected FIRFilterImplementation firFilter = new FIRFilterImplementation(3);
    protected IIRFilterImplementation iirFilter = new IIRFilterImplementation(2);

    protected double[] a = new double[3];
    protected double[] b = new double[2];
    protected double gain;

    public HighpassFilterButterworthSection(double cutoffFrequencyHz, double k, double n, double Fs)
    {
        // pre-warp omegac and invert it
        double omegac = 1.0 /(2.0 * Fs * Math.Tan(Math.PI * cutoffFrequencyHz / Fs));

        // compute zeta
        double zeta = -Math.Cos(Math.PI * (2.0 * k + n - 1.0) / (2.0 * n));

        // fir section
        this.a[0] = 4.0 * Fs * Fs;
        this.a[1] = -8.0 * Fs * Fs;
        this.a[2] = 4.0 * Fs * Fs;

        //iir section
        //normalize coefficients so that b0 = 1
        //and higher-order coefficients are scaled and negated
        double b0 = (4.0 * Fs * Fs) + (4.0 * Fs * zeta / omegac) + (1.0 / (omegac * omegac));
        this.b[0] = ((2.0 / (omegac * omegac)) - (8.0 * Fs * Fs)) / (-b0);
        this.b[1] = ((4.0 * Fs * Fs)
                   - (4.0 * Fs * zeta / omegac) + (1.0 / (omegac * omegac))) / (-b0);
        this.gain = 1.0 / b0;
    }

    public double compute(double input)
    {
        // compute the result as the cascade of the fir and iir filters
        return this.iirFilter.compute
            (this.firFilter.compute(this.gain * input, this.a), this.b);
    }
}

using System;

namespace DSP
{
    public class BandpassFilterButterworthImplementation
    {
        protected LowpassFilterButterworthImplementation lowpassFilter;
        protected HighpassFilterButterworthImplementation highpassFilter;
        
        public BandpassFilterButterworthImplementation
           (double bottomFrequencyHz, double topFrequencyHz, int numSections, double Fs)
        {
            this.lowpassFilter = new LowpassFilterButterworthImplementation
                                 (topFrequencyHz, numSections, Fs);
            this.highpassFilter = new HighpassFilterButterworthImplementation
                                  (bottomFrequencyHz, numSections, Fs);
        }

        public double compute(double input)
        {
            // compute the result as the cascade of the highpass and lowpass filters
            return this.highpassFilter.compute(this.lowpassFilter.compute(input));
        }
    }
}

using System;

namespace DSP
{
    public class FIRFilterImplementation
    {
        protected double[] z;
        public FIRFilterImplementation(int order)
        {
            this.z = new double[order];
        }

        public double compute(double input, double[] a)
        {
            // computes y(t) = a0*x(t) + a1*x(t-1) + a2*x(t-2) + ... an*x(t-n)
            double result = 0;

            for (int t = a.Length - 1; t >= 0; t--)
            {
                if (t > 0)
                {
                    this.z[t] = this.z[t - 1];
                }
                else
                {
                    this.z[t] = input;
                }
                result += a[t] * this.z[t];
            }
            return result;
        }
    }
}

namespace DSP
{
    public class IIRFilterImplementation
    {
        protected double[] z;
        public IIRFilterImplementation(int order)
        {
            this.z = new double[order];
        }

        public double compute(double input, double[] a)
        {
            // computes y(t) = x(t) + a1*y(t-1) + a2*y(t-2) + ... an*y(t-n)
            // z-transform: H(z) = 1 / (1 - sum(1 to n) [an * y(t-n)])
            // a0 is assumed to be 1
            // y(t) is not stored, so y(t-1) is stored at z[0], 
            // and a1 is stored as coefficient[0]

            double result = input;

            for (int t = 0; t < a.Length; t++)
            {
                result += a[t] * this.z[t];
            }
            for (int t = a.Length - 1; t >= 0; t--)
            {
                if (t > 0)
                {
                    this.z[t] = this.z[t - 1];
                }
                else
                {
                    this.z[t] = result;
                }
            }
            return result;
        }
    }
}

History

• 9^th June, 2019: Initial submission