计算2阶Butterworth低通滤波器系数

Question

有了以下参数，您该如何计算以下差分方程的系数？

采样率：10kHz
截止频率：1kHz

我知道差分方程将采用以下形式，但不知道如何实际计算并得到系数b0、b1、b2、a1、a2的数字:

y(n)  =  b0.x(n) + b1.x(n-1) + b2.x(n-2) + a1.y(n-1) + a2.y(n-2)

我最终将在C++中实现这个低通滤波器，但我需要知道如何计算系数，然后才能开始任何工作。

- Daniel

---

这个问题似乎与编程无关，而是涉及信号处理理论。 - talonmies

所以，问一个显而易见的问题，你有查看相关的维基百科条目吗？ - Karoly Horvath

您提供的公式看起来像是一个通用的二阶微分方程，但您并没有提供边界条件（或等价条件），因此目前这个问题似乎无法回答。您能否提供更多上下文信息？ - Dan Nissenbaum

@talonmies 我以前用过 Stack Overflow 并且得到了很好的回复，但我不确定其他哪里可以问。 - Daniel

@KarolyHorvath 是的，我已经查过维基百科了，那是我第一个去的地方。 - Daniel

Answer 1

给你。ff是频率比，在你的情况下为0.1：

    const double ita =1.0/ tan(M_PI*ff);
    const double q=sqrt(2.0);
    b0 = 1.0 / (1.0 + q*ita + ita*ita);
    b1= 2*b0;
    b2= b0;
    a1 = 2.0 * (ita*ita - 1.0) * b0;
    a2 = -(1.0 - q*ita + ita*ita) * b0;

结果为：

b0=0.0674553
b1=0.134911
b2=0.0674553
a1=1.14298
a2=-0.412802

- pentadecagon

---

这有帮助，我能否修改这个程序来计算不同的截止频率和采样频率，以获得任何输入的系数？ - Daniel

只需使用提供的公式，其中ff是截止频率和采样频率之比：ff = f_cutoff / f_sampling。 - pentadecagon

这里有更多关于计算巴特沃斯滤波器增益和编码巴特沃斯滤波器的内容：http://www.centerspace.net/blog/nmath/iir-filtering-with-butterworth-filters/ - Paul

我没有信号/电子工程背景。请问有人可以解释一下为什么系数不必加起来等于1吗？看起来OP发布的方程式会输入x[n]、x[n-1]等，然后输出y[n]，中心路径与一个完全不同的“平均”信号路径相对应。 - Branden Keck

这些数字加起来确实等于1，b0+b1+b2+a1+a2=1，所以我不确定你的意思是什么。 - pentadecagon

Answer 2

对于那些想知道其他答案中的神奇公式来自哪里的人，这里是按照这个例子推导的过程。

从Butterworth滤波器的传递函数开始：

G(s) = wc^2 / (s^2 + s*sqrt(2)*wc + wc^2)

其中，wc是截止频率，应用双线性z变换，即将s = 2/T*(1-z^-1)/(1+z^-1)代入：

G(z) = wc^2 / ((2/T*(1-z^-1)/(1+z^-1))^2 + (2/T*(1-z^-1)/(1+z^-1))*sqrt(2)*wc + wc^2)

T是采样周期[s]。

截止频率需要进行预弯曲以补偿由z变换引入的模拟和数字频率之间的非线性关系：

wc = 2/T * tan(wd*T/2)

其中wd是所需的截止频率[rad/s]。

为了方便起见，令C = tan(wd*T/2)，使得wc = 2/T*C。

将它代入方程，2/T因子消失：

G(z) = C^2 / ((1-z^-1)/(1+z^-1))^2 + (1-z^-1)/(1+z^-1)*sqrt(2)*C + C^2)

将分子和分母乘以(1+z^-1)^2并展开，得到：

G(z) = C^2*(1 + 2*z^-1 + z^-2) / (1 + sqrt(2)*C + C^2 + 2*(C^2-1)*z^-1 + (1-sqrt(2)*C+C^2)*z^-2')

现在，将分子和分母都除以分母的常数项。为方便起见，令D = 1 + sqrt(2)*C + C^2：

G(z) = C^2/D*(1 + 2*z^-1 + z^-2) / (1 + 2*(C^2-1)/D*z^-1 + (1-sqrt(2)*C+C^2)/D*z^-2')

这个形式等价于我们要找的形式：

G(z) = (b0 + b1*z^-1 + b2*z^-1) / (1 + a1*z^-1 +a2*z^-2)

因此，我们通过等式得到系数：

a0 = 1

a1 = 2*(C^2-1)/D

a2 = (1-sqrt(2)*C+C^2)/D

b0 = C^2/D

b1 = 2*b0

b2 = b0

上述公式中，D = 1 + sqrt(2)*C + C^2，C = tan(wd*T/2)，其中wd为所需截止频率（单位：弧度/秒），T为采样周期（单位：秒）。请注意保留HTML标签。

- djvg

Answer 3

如果你需要一个高通滤波器系数，你只需要使用相同的代码：

const double ita =1.0/ tan(M_PI*ff);
const double q=sqrt(2.0);
b0 = 1.0 / (1.0 + q*ita + ita*ita);
b1= 2*b0;
b2= b0;
a1 = 2.0 * (ita*ita - 1.0) * b0;
a2 = -(1.0 - q*ita + ita*ita) * b0;

但是在将所有的b项乘以ita^2并取反b1之后。

b0 = b0*ita*ita;
b1 = -b1*ita*ita;
b2 = b2*ita*ita;

现在您拥有了一个二阶高通滤波器

- wmohr

---

澄清：a1和a2保持不变，只有b的值会改变。 - velochy

Answer 4

您可以使用这个链接获取特定采样率和截止频率的n阶Butterworth滤波器系数。为了测试结果，您可以使用MATLAB获取系数并与程序的输出进行比较。 http://www.exstrom.com/journal/sigproc

fnorm = f_cutoff/(f_sample_rate/2); % normalized cut off freq, http://www.exstrom.com/journal/sigproc
% Low pass Butterworth filter of order N
[b1, a1] = butter(nth_order, fnorm,'low');

- Richard Le

---

请注意，Matlab使用的Fnorm和SCIPY使用的Wn的值是Pentadecagon答案中使用的FF值的一半。我浪费了半天时间没有注意到这一点... - Tunneller

Answer 5

这里我的a和b被搞混了，但是我的代码看起来像这样:

//Boulanger and Lazzarini, The Audio Programming Book, pg 484
void calculateDifferenceEquation() {
    float lambda = 1.0 / (tan(M_PI * mFrequency / mSampleRate));
    a0 = 1.0 / (1.0 + (2.0 * lambda) + pow(lambda, 2.0));
    a1 = 2.0 * a0;
    a2 = a0;
    b1 = 2.0 * a0 * (1.0 - pow(lambda, 2.0));
    b2 = a0 * (1.0 - (2.0 * lambda) + pow(lambda, 2.0));
}

- Clay Bailey

---

嗨，你的a和b被交换了，但是与下面pentadecagon发布的相比，你的b带有负号。实际上，似乎存在一个我无法解决的微妙差别。在下面的数学中，sqrt(2)lambda似乎是正确的，而你的方程式中是2lambda。我输入了你的方程式，它确实平滑了，但对于相同的lambda值，它更加激进。 - Tunneller

Answer 6

最好的方法是使用类似LabVIEW的工具来模拟您的滤波器，根据您的fc和fs获取系数，然后在C语言中使用这些系数。最后将代码烧入您的微控制器，并将响应与LabVIEW或Simulink中的响应进行比较。

- Bharat

---

嘿，Bharat，感谢您的回复。您能否编辑帖子并添加一些示例？ - StackExchange What The Heck