Click here to Skip to main content
16,020,261 members
Articles / General Programming / Algorithms
Alternative
Tip/Trick

Lowpass, Highpass, and Bandpass Butterworth Filters in C#

Rate me:
Please Sign up or sign in to vote.
4.45/5 (8 votes)
9 Jun 2019CPOL1 min read 42.5K   2.3K   13   7
Implementations of lowpass, highpass, and bandpass nth-order Butterworth filters in C#, with design documentation.

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 nth-order FIR filter for the zero (numerator) polynomials and an implementation of an nth-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.

C#
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);
    }
}

 

C#
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);
    }
}

 

C#
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));
        }
    }
}
C#
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;
        }
    }
}
C#
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

  • 9th June, 2019: Initial submission

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)



Comments and Discussions

 
QuestionFs is a rather confusing parameter Pin
Member 1014767325-May-23 5:37
Member 1014767325-May-23 5:37 
QuestionFYI - The download is missing the actual source code Pin
Lee.W.Spencer25-Feb-22 2:44
Lee.W.Spencer25-Feb-22 2:44 
QuestionQuick and Dirty Frontend Pin
Pookalton19-Jan-21 20:47
Pookalton19-Jan-21 20:47 
QuestionA Suggested Method To Add Pin
BKrell23-Nov-20 23:30
BKrell23-Nov-20 23:30 
QuestionHow can I use this on a real signal? Please add explanation. Pin
Hamed Nik12-Sep-19 0:02
Hamed Nik12-Sep-19 0:02 
AnswerRe: How can I use this on a real signal? Please add explanation. Pin
HardworkingPearl5-Jan-20 21:10
HardworkingPearl5-Jan-20 21:10 
SuggestionPlease add a Winform and a couple of nice signal sources for a real demo Pin
asiwel11-Jun-19 12:04
professionalasiwel11-Jun-19 12:04 
Hi, Jim. This looks like a very nice article and toolbox resource. And, of course, I went back and re-read the Wikipedia article on Butterworth filters. And your math document was a nice touch.
Butterworth filter - Wikipedia[^]

But for myself and other readers, I'd like to see a real demo program; for instance, a simple Winform, a couple of typical signal sources, and a few pretty graphs or charts similar to the ones shown on the Wiki page. Then I could have fun playing with these filters, changing the signal source characteristics, adding a few unit tests, etc.

In all likelihood, you already have done this sort of thing during the process of developing and testing your classes, so this would seem to be a simple request. Thank you.

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.