C# 中的低通、高通和带通巴特沃斯滤波器
4.45/5 (8投票s)
C# 中任意阶数的低通、高通和带通巴特沃斯滤波器的实现，附带设计文档。
引言
这是一个 C# 实现的数字低通、高通和带通巴特沃斯滤波器，阶数任意（n 个级联的二阶部分）。
背景
包含一份 Word 文档，介绍了通过双线性 z 变换进行滤波器设计。
设计从一个连续时间低通巴特沃斯滤波器开始，并使用双线性 z 变换推导出等效数字滤波器的系数。然后，变量替换 s=1/s 将低通滤波器转换为高通滤波器，并使用双线性 z 变换推导出等效的高通数字滤波器的系数。
在代码中，低通和高通滤波器是根据设计文档中的系数推导来实现的。带通滤波器被实现为低通滤波器和高通滤波器的级联。
代码包括一个 n 阶 FIR 滤波器，用于零（分子）多项式，以及一个 n 阶 IIR 滤波器，用于极点（分母）多项式。
Using the Code
代码分为五个文件
- LowpassFilterButterworthImplementation.cs
- HighpassFilterButterworthImplementation.cs
- BandpassFilterButterworthImplementation.cs
- FIRFilterImplementation.cs
- IIRFilterImplementation.cs
要使用 LowpassFilterButterworthImplementation、HighpassFilterButterworthImplementation 或 BandpassFilterButterworthImplementation，请使用所需的截止频率或频率、二阶部分数量和采样频率 (Fs) 构造对象。然后，反复调用对象的 compute() 方法，传入一个输入样本并接收一个过滤后的样本作为返回值。
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;
}
}
}
历史
- 2019 年 6 月 9 日：初始提交