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:
- LowpassFilterButterworthImplementation.cs
- HighpassFilterButterworthImplementation.cs
- BandpassFilterButterworthImplementation.cs
- FIRFilterImplementation.cs
- 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)
{
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));
this.a[0] = omegac * omegac;
this.a[1] = 2.0 * omegac * omegac;
this.a[2] = omegac * omegac;
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)
{
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)
{
double omegac = 1.0 /(2.0 * Fs * Math.Tan(Math.PI * cutoffFrequencyHz / Fs));
double zeta = -Math.Cos(Math.PI * (2.0 * k + n - 1.0) / (2.0 * n));
this.a[0] = 4.0 * Fs * Fs;
this.a[1] = -8.0 * Fs * Fs;
this.a[2] = 4.0 * Fs * Fs;
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)
{
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)
{
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)
{
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)
{
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