在C++中实现二阶低通滤波器，如何计算系数？

Question

我正在努力寻找生成低通滤波器系数的合适算法。 我写了以下代码，使用另一个 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]

这又不是很合适。

- nowox

---

如果这与C++有任何特定关系，那么计算浮点值时使用的精度可能会有所不同：每个程序员都应该了解的浮点运算知识 - πάντα ῥεῖ

你可以查看 Embree 和 Danieli 编写的 C++ 数字信号处理算法书籍。涵盖了 IIR 和 FIR（如果需要线性相位 FIR，则是最佳选择）。该书附带代码。 - doug

Answer 1

如果z变换后的传递函数被定义为：

H(z) = ( b0 + b1/z + b2/z^2 ) / ( a0 + a1/z + a2/z^2 )

那么我认为系数a1和a2的公式产生了错误符号的值。（至少，当我重新计算时发现了这一点，但我可能是错的）。

以下代码给出了a1和a2相反的符号。

#include <iostream>
#include <iomanip>
#include <complex>
#include <vector>
#include <cmath>
using namespace std;

const double PI = 4.0 * atan( 1.0 );
const complex<double> iPI( 0.0, PI );
const double SQRT2 = sqrt( 2.0 );

// const double factor = 0.5;    // use for Matlab comparison
const double factor = 1.0;       // as original post

//----------------------------------------------------------------------

void getCoefficients( double fd, double fs, vector<double> &a, vector<double> &b )
{
   double alpha = 1.0 / tan( factor * PI * fd / fs );
   b[0] = 1.0 / ( 1.0 + SQRT2 * alpha + alpha * alpha );
   b[1] = 2 * b[0];
   b[2] = b[0];
   a[0] = 1.0;
   a[1] = 2.0 * ( 1.0 - alpha * alpha ) * b[0];
   a[2] = ( 1.0 - SQRT2 * alpha + alpha * alpha ) * b[0];
}

//----------------------------------------------------------------------

void bode( const vector<double> &a, const vector<double> &b, double f, double fs, double &mag, double &phase )
{
   complex<double> invz = ( 1.0 - iPI * factor * f / fs ) / ( 1.0 + iPI * factor * f / fs );
   complex<double> H = ( b[0] + b[1] * invz + b[2] * invz * invz ) / ( a[0] + a[1] * invz + a[2] * invz * invz );
   mag = abs( H );
   phase = atan2( H.imag(), H.real() );
}

//----------------------------------------------------------------------

int main()
{
   vector<double> a(3), b(3);
   double fs = 10000;
   double fd = 1000;

   // Get Butterworth coefficients in H(z) = ( b0 + b1/z + b2/z^2 ) / ( a0 + a1/z + a2/z^2 )
   getCoefficients( fd, fs, a, b );
   cout << "a:\t";
   for ( double e : a ) cout << e << '\t';
   cout << "\nb:\t";
   for ( double e : b ) cout << e << '\t';
   cout << "\n\n";

   double logmin = 1.0, logmax = 4.0;
   int nw = 31;
   double dlog = ( logmax - logmin ) / ( nw - 1 );
   #define FMT << scientific << setw(15) << setprecision(4) <<

   cout FMT "f(Hz)" FMT "mag" FMT "db" FMT "phase(deg)" << '\n';
   for ( int i = 0; i < nw; i++ )
   {
      double f, mag, phase;
      f = pow( 10.0, logmin + i * dlog );
      bode( a, b, f, fs, mag, phase );

      double phasedeg = phase * 180.0 / PI;
      double db = 20.0 * log( mag ) / log( 10.0 );
      cout FMT f FMT mag FMT db FMT phasedeg << '\n';
   }
}

输出：

a:  1   -1.14298    0.412802    
b:  0.0674553   0.134911    0.0674553   

          f(Hz)            mag             db     phase(deg)
     1.0000e+01     1.0000e+00    -3.7956e-08    -7.8347e-01
     1.2589e+01     1.0000e+00    -9.5341e-08    -9.8635e-01
     1.5849e+01     1.0000e+00    -2.3949e-07    -1.2418e+00
     1.9953e+01     1.0000e+00    -6.0156e-07    -1.5634e+00
     2.5119e+01     1.0000e+00    -1.5111e-06    -1.9683e+00
     3.1623e+01     1.0000e+00    -3.7956e-06    -2.4783e+00
     3.9811e+01     1.0000e+00    -9.5341e-06    -3.1205e+00
     5.0119e+01     1.0000e+00    -2.3949e-05    -3.9296e+00
     6.3096e+01     9.9999e-01    -6.0156e-05    -4.9494e+00
     7.9433e+01     9.9998e-01    -1.5110e-04    -6.2354e+00
     1.0000e+02     9.9996e-01    -3.7954e-04    -7.8588e+00
     1.2589e+02     9.9989e-01    -9.5331e-04    -9.9113e+00
     1.5849e+02     9.9972e-01    -2.3942e-03    -1.2513e+01
     1.9953e+02     9.9931e-01    -6.0114e-03    -1.5821e+01
     2.5119e+02     9.9826e-01    -1.5084e-02    -2.0052e+01
     3.1623e+02     9.9566e-01    -3.7791e-02    -2.5501e+01
     3.9811e+02     9.8920e-01    -9.4310e-02    -3.2581e+01
     5.0119e+02     9.7352e-01    -2.3312e-01    -4.1849e+01
     6.3096e+02     9.3720e-01    -5.6339e-01    -5.3957e+01
     7.9433e+02     8.6132e-01    -1.2967e+00    -6.9313e+01
     1.0000e+03     7.3050e-01    -2.7276e+00    -8.7273e+01
     1.2589e+03     5.5943e-01    -5.0451e+00    -1.0563e+02
     1.5849e+03     3.9180e-01    -8.1387e+00    -1.2189e+02
     1.9953e+03     2.5949e-01    -1.1718e+01    -1.3493e+02
     2.5119e+03     1.6715e-01    -1.5538e+01    -1.4496e+02
     3.1623e+03     1.0636e-01    -1.9464e+01    -1.5262e+02
     3.9811e+03     6.7339e-02    -2.3435e+01    -1.5850e+02
     5.0119e+03     4.2546e-02    -2.7423e+01    -1.6305e+02
     6.3096e+03     2.6859e-02    -3.1418e+01    -1.6660e+02
     7.9433e+03     1.6951e-02    -3.5416e+01    -1.6939e+02
     1.0000e+04     1.0696e-02    -3.9415e+01    -1.7159e+02

关于Matlab...如果你将代码顶部的因子更改为注释掉的值（0.5），那么你将得到Matlab的值：

a:      1       -1.56102        0.641352        
b:      0.0200834       0.0401667       0.0200834  

 

- lastchance

Answer 2

我没有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帖子相比，你的"a"值具有负号，我没有检查这对结果的影响。

- Tunneller

Answer 3

对于Butterworth，您可以使用以下内容：

void butter(float fc, float fs, bool lowPass = true) {
    const float k = tanf(M_PI * fc/(fs*2));
    const float k2 = k * k;
    const float q = 1 / sqrtf(2);
    const float norm = 1 / (1 + k / q + k2);
    if (lowPass) {
        b0 = k2 * norm;
        b1 = 2 * b0;
        b2 = b0;
    } else {
        b0 = 1 * norm;
        b1 = -2 * b0;
        b2 = b0;
    }
    a1 = 2 * (k2 - 1) * norm;
    a2 = (1 - k / q + k2) * norm;
}

在Matlab中，你可以使用来获得完全相同的系数。

- nowox