计算二阶巴特沃斯低通滤波器的系数答案

【问题标题】：Calculate Coefficients of 2nd Order Butterworth Low Pass Filter计算二阶巴特沃斯低通滤波器的系数
【发布时间】：2014-01-22 09:00:35
【问题描述】：

随着，

采样频率：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++ 中实现这个 LPF，但我需要先知道如何实际计算系数，然后才能使用它

【问题讨论】：

这个问题似乎是题外话，因为它是关于信号处理理论而不是编程
那么，很明显，您检查过相关的维基百科条目吗？
您提供的公式看起来像一个通用的二阶微分方程，并且您没有提供边界条件（或等效条件），因此目前这个问题似乎无法回答。你能提供更多背景信息吗？
@talonmies 我之前使用过堆栈溢出并收到了很好的响应，但我不确定我还能问到哪里
@KarolyHorvath 是的，我查过维基百科，我去的第一个地方

标签： c++ signal-processing

【解决方案1】：

我的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));
}

【讨论】：

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

【解决方案2】：

对于那些想知道其他答案的神奇公式来自哪里的人，这里是this example 之后的推导。

从巴特沃斯滤波器的传递函数开始

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 是所需的截止频率 [rad/s]，T 是采样周期 [s]。

【讨论】：

【解决方案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;

现在你有一个二阶高通滤波器

【讨论】：

澄清：a1 和 a2 保持不变。只有 b 值改变

【解决方案4】：

最好的方法是使用诸如实验室视图之类的东西来模拟您的滤波器并根据您的 fc 和 fs 获取系数。然后在 c 中使用它们。最后将代码烧录到您的 microcon 中。并将响应与实验室视图或 simulink 中的响应进行比较。

【讨论】：

嘿巴拉特，感谢您的回复。您可以编辑帖子并添加一些示例吗？

【解决方案5】：

您可以使用此链接获取具有特定采样率和频率截止的 n 阶巴特沃斯滤波器的系数。为了测试结果。您可以使用 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');

【讨论】：

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

【解决方案6】：

给你。 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

【讨论】：

这有帮助，有没有办法可以修改它以便能够计算不同的截止和采样频率以获得任何输入的系数？
使用提供的公式即可，ff是截止和采样频率的比值：ff=f_cutoff / f_sampling
这是一个bit more，关于计算巴特沃斯滤波器增益和编码巴特沃斯滤波器。
我不是信号/EE 背景。有人可以解释为什么系数不必加到 1 吗？看起来 OP 发布的方程式将采用 x[n]、x[n-1] 等，并以完全不同的“平均”信号路径为中心输出 y[n]
数字加起来确实是 1，b0+b1+b2+a1+a2=1，所以我不确定你的意思。