计算二阶巴特沃斯低通滤波器的系数
Calculate Coefficients of 2nd Order Butterworth Low Pass Filter
使用,
采样频率: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,但我需要知道如何首先实际计算系数,然后才能使用它
开始吧。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.144298
a2=-0.412802
对于那些想知道其他答案中的神奇公式来自哪里的人,下面是这个例子的推导。
从巴特沃斯滤波器的传递函数开始
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]。
您可以使用此链接获得具有特定采样率和频率切割的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');
FYI如果你需要一个高通滤波器系数,你所要做的就是使用相同的代码:
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;
现在你有一个二阶高通滤波器
我的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));
}
最好的方法是使用类似实验室视图的东西来模拟滤波器,并根据fc和fs获得系数。然后在c中使用它们,最后把代码烧到你的微码中。并将响应与实验室视图或simulink中的响应进行比较。
- 为什么"do while"循环不断退出,即使条件计算结果为 false?
- 表示"accepting anything for this template argument" C++概念的通配符
- 递归函数计算序列中的平方和(并输出过程)
- (C++)分析树以计算返回错误值的简单算术表达式
- 我的字符计数代码计算错误.为什么
- 在计算中使用二的幂有多有利可图
- 如何计算该程序的复杂性?是否有任何其他复杂性较低的解决方案
- 计算二阶巴特沃斯低通滤波器的系数
- 对于给定的n个数字序列,找到具有最大总和和尽可能低的计算复杂度的子串
- 计算iir滤波器的w个系数
- 带通FIR滤波器
- C++中的带通巴特沃兹滤波器实现
- 带FFT卷积的低通FIR滤波器-重叠加法,原因和方式
- 计算梅尔滤波器组系数的公式
- 如何使用C++计算海尔斯通序列中的奇数个数
- 以低复杂度计算K个连续数字模组1000000007的LCM的程序
- C++,计算具有多个"&&"且没有较低优先级运算符的表达式
- 用FFTW在C语言中实现FFT低通滤波器
- 计算与通配符匹配的下一个最近日期时间
- IIR低通滤波器在c++中的实现