我正在努力寻找合适的算法来生成低通滤波器的系数。我从another SO问题中提取了butterworthLowPass
代码,并编写了以下代码:
class Filter {
float fs;
float a1, a2, b0, b1, b2;
float v1, v2;
public:
Filter(float fs) : fs(fs) {}
float compute(float x)
{
float v0 = x - a1 * v1 - a2 * v2;
float y = b0 * v0 + b1 * v1 + b2 * v2;
v2 = v1;
v1 = v0;
return y;
}
void butterworthLowPass(float fc)
{
const float ita = 1.0 / tan(M_PI * fc / fs);
const float 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;
}
void chebyshevLowPass(float f, float ripple)
{
// TODO: implement
}
void ellipticLowPass(float f, float ripple, float stopband)
{
// TODO: implement
}
void displayCoefficients()
{
fprintf(stderr, "a = [ %f, %f, %f ]\n", 1.f, a1, a2);
fprintf(stderr, "b = [ %f, %f, %f ]\n", b0, b1, b2);
}
};
这里的主要:
int main()
{
const float fs = 10000;
Filter biquad(fs);
biquad.butterworthLowPass(1000);
biquad.displayCoefficients();
}
当计算系数并使用Matlab显示滤波器响应时,我得到了-15dB增益,这不是我所期望的。
a = [ 1.000000, 1.142980, -0.412802 ];
b = [ 0.067455, 0.134911, 0.067455 ];
[mag, phase, wout] = bode(tf(b,a,1/(fs/2)));
subplot(2,1,1);
semilogx(wout(:,1)/(2*pi), 20*log10(squeeze(mag)), '-b'); zoom on; grid on;
axis tight
ylim([-40 0])
title('magnitude'); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
subplot(2,1,2);
semilogx(wout(:,1)/(2*pi), squeeze(phase), '-r'); zoom on; grid on;
axis tight
Matlab给了我一些不同的东西:
>> fs = 10000;
>> fc = 1000;
>> [b, a] = butter(2, 1000/10000)
b =
0.0201 0.0402 0.0201
a =
1.0000 -1.5610 0.6414
进一步实验
我试着克隆ruohoruotsi/Butterworth-Filter-Design.git
并测试:
#include "Butterworth.h"
using namespace std;
int main() {
float fs = 20000;
float fc = 500;
int order = 2;
vector<Biquad> coeffs; // second-order sections (sos)
Butterworth butterworth;
bool designedCorrectly = butterworth.loPass(fs, fc, 0, order, coeffs, 1.0f);
std::cout << "Designed correctly? " << designedCorrectly << std::endl;
std::cout << "a = [" << 1.f << ", " << coeffs[0].a1 << ", " << coeffs[0].a2
<< "]" << std::endl;
std::cout << "b = [" << coeffs[0].b0 << ", " << coeffs[0].b1 << ", "
<< coeffs[0].b2 << "]" << std::endl;
}
它给了我:
a = [1, 1.77863, -0.800803]
b = [1, 2, 1]
这又不太合适。
3条答案
按热度按时间cngwdvgl1#
如果z变换后的传递函数定义为
那么我认为系数a1和a2的公式产生了错误符号的值。(至少,这是我在重新工作时发现的,但我可能错了)。
下面的代码给出了a1和a2的相反符号。
输出:
对于Matlab…好吧,如果你把代码顶部的因子改为注解掉的值(0.5),那么你会得到Matlab的值:
a7qyws3x2#
我没有matlab,但是python。
butter(2,.2) = [0.06745527, 0.13491055, 0.06745527], [ 1., -1.1429805, 0.4128016]
和
butter(2,.1) = [0.02008337, 0.04016673, 0.02008337],[ 1. , -1.56101808, 0.64135154]
所以你的Matlab值在FF上偏离了“2”的因子(我曾经犯过同样的错误,并在Calculate Coefficients of 2nd Order Butterworth Low Pass Filter上发表了评论)
“Butterworth.h”代码通过乘以
(1.0 + q*ita + ita*ita)
对“b”值进行了归一化。与Matlab、Scipy和原始SO post使用的值相比,您的“a”值有负号,我也没有检查它的贡献。omvjsjqw3#
对于巴特沃思,您可以使用以下命令:
在Matlab中,您可以使用
butter(2, fc/fs)
获得完全相同的系数