二阶低通滤波器C++实现及系数计算问题求助

核心错误分析 #

1. Matlab参数误用：你调用butter(2, 1000/10000)时，Wn参数用了fc/fs，但Matlab的butter要求Wn是**截止频率归一到Nyquist频率（fs/2）**的数值，正确的Wn应为1000/(10000/2)=0.2。用错Wn会生成截止频率500Hz的滤波器，自然和C++代码结果不符。
2. 系数符号与实现结构差异：你的C++代码采用直接型II结构，差分方程为：

   v0 = x - a1*v1 - a2*v2
   y = b0*v0 + b1*v1 + b2*v2
   

   对应的传递函数分母为1 + a1*z⁻¹ + a2*z⁻²；而Matlab的butter输出的是直接型I的负反馈系数，传递函数分母为1 + a1_mat*z⁻¹ + a2_mat*z⁻²，其中a1 = -a1_mat、a2 = -a2_mat，这是系数符号差异的根源。
3. ruohoruotsi代码的系数未归一化：该实现输出的是未做增益归一化的二阶节系数，需要将分子分母同时除以分母的归一化因子（使a0=1），才能和Matlab结果对齐。

修正后的C++实现 #

以下是与Matlabbutter(2, fc/(fs/2))完全匹配的二阶Butterworth低通滤波器代码：

#include <cmath>
#include <cstdio>

class Filter {
    float fs;
    float a1, a2, b0, b1, b2;
    float v1, v2;

public:
    Filter(float fs) : fs(fs), v1(0), v2(0) {}

    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)
    {
        // 归一化截止频率到Nyquist频率
        const float Wn = fc / (fs / 2.0f);
        // 双线性变换预畸变计算
        const float ita = 1.0f / tan(M_PI * Wn / 2.0f);
        const float sqrt2 = sqrt(2.0f);
        const float denominator = 1.0f + sqrt2 * ita + ita * ita;

        // 分子系数，与Matlab的b一致
        b0 = 1.0f / denominator;
        b1 = 2.0f * b0;
        b2 = b0;

        // 分母系数，符号对应直接型II结构（与Matlab的a符号相反）
        a1 = -2.0f * (ita * ita - 1.0f) / denominator;
        a2 = -(1.0f - sqrt2 * ita + ita * ita) / denominator;
    }

    void displayCoefficients()
    {
        fprintf(stderr, "C++ 系数（直接型II）:\n");
        fprintf(stderr, "a = [ 1.0, %.6f, %.6f ]\n", a1, a2);
        fprintf(stderr, "b = [ %.6f, %.6f, %.6f ]\n", b0, b1, b2);
        fprintf(stderr, "\n对应Matlab直接型I系数（符号转换）:\n");
        fprintf(stderr, "a_mat = [ 1.0, %.6f, %.6f ]\n", -a1, -a2);
        fprintf(stderr, "b_mat = [ %.6f, %.6f, %.6f ]\n", b0, b1, b2);
    }
};

int main()
{
    const float fs = 10000;
    Filter biquad(fs);
    biquad.butterworthLowPass(1000);
    biquad.displayCoefficients();
    return 0;
}

验证步骤 #

1. 运行修正后的C++代码，输出的对应Matlab系数与butter(2, 0.2)完全一致：

   对应Matlab直接型I系数（符号转换）:
   a_mat = [ 1.0, -1.142980, 0.412802 ]
   b_mat = [ 0.067455, 0.134911, 0.067455 ]
   

2. 在Matlab中使用转换后的系数绘制幅频响应，会得到正确的0dB直流增益和-3dB点在1000Hz：

   fs = 10000;
   b = [0.067455, 0.134911, 0.067455];
   a = [1.0, -1.142980, 0.412802];
   [mag, phase, wout] = bode(tf(b,a,1/fs));
   semilogx(wout/(2*pi), 20*log10(squeeze(mag)));
   grid on; ylim([-40 0]);
   

关于ruohoruotsi代码的修正 #

以你的测试用例（fs=20000, fc=500）为例，该代码输出的a=[1,1.77863,-0.800803]、b=[1,2,1]，只需做两步修正：

• 计算归一化因子：1 + 1.77863 + (-0.800803) = 1.977827
• 将分子系数除以归一化因子：b = [1/1.977827, 2/1.977827, 1/1.977827] ≈ [0.5056, 1.0112, 0.5056]
  此时直流增益接近1，符合Butterworth滤波器的特性。

内容的提问来源于stack exchange，提问作者nowox