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

huangapple
go评论99阅读模式
英文:

Implement 2nd order low pass filter in C++, how to compute coefficients?

问题

I see that you're working on designing and implementing various low-pass filter algorithms in C++. How can I assist you further with this code or any specific questions you have about it?

英文:

I am struggling to find the proper algorithms for generating the coefficients for low pass filters. I wrote the following taking the butterworthLowPass code from another SO question:

  1. class Filter {
  2.     float fs;
  3.     float a1, a2, b0, b1, b2;
  4.     float v1, v2;
  6.    public:
  7.     Filter(float fs) : fs(fs) {}
  9.     float compute(float x)
  10.     {
  11.         float v0 = x - a1 * v1 - a2 * v2;
  12.         float y = b0 * v0 + b1 * v1 + b2 * v2;
  13.         v2 = v1;
  14.         v1 = v0;
  15.         return y;
  16.     }
  18.     void butterworthLowPass(float fc)
  19.     {
  20.         const float ita = 1.0 / tan(M_PI * fc / fs);
  21.         const float q = sqrt(2.0);
  22.         b0 = 1.0 / (1.0 + q * ita + ita * ita);
  23.         b1 = 2 * b0;
  24.         b2 = b0;
  25.         a1 = 2.0 * (ita * ita - 1.0) * b0;
  26.         a2 = -(1.0 - q * ita + ita * ita) * b0;
  27.     }
  29.     void chebyshevLowPass(float f, float ripple)
  30.     {
  31.         // TODO: implement
  32.     }
  34.     void ellipticLowPass(float f, float ripple, float stopband)
  35.     {
  36.         // TODO: implement
  37.     }
  39.     void displayCoefficients()
  40.     {
  41.         fprintf(stderr, "a = [ %f, %f, %f ]\n", 1.f, a1, a2);
  42.         fprintf(stderr, "b = [ %f, %f, %f ]\n", b0, b1, b2);
  43.     }
  44. };

Here the main:

  1. int main()
  2. {
  3.     const float fs = 10000;
  4.     Filter biquad(fs);
  5.     biquad.butterworthLowPass(1000);
  6.     biquad.displayCoefficients();
  7. }

When taking the computed coefficients and using Matlab to display the filter response, I get a -15dB gain which is not what I expect.

  1. a = [ 1.000000, 1.142980, -0.412802 ];
  2. b = [ 0.067455, 0.134911, 0.067455 ];
  3. [mag, phase, wout] = bode(tf(b,a,1/(fs/2)));
  4. subplot(2,1,1);
  5. semilogx(wout(:,1)/(2*pi), 20*log10(squeeze(mag)), '-b'); zoom on; grid on;
  6. axis tight
  7. ylim([-40 0])
  8. title('magnitude'); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
  9. subplot(2,1,2);
  10. semilogx(wout(:,1)/(2*pi), squeeze(phase), '-r'); zoom on; grid on;
  11. axis tight

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

Matlab gives me something very different:

  1. >> fs = 10000;
  2. >> fc = 1000;    
  3. >> [b, a] = butter(2, 1000/10000)
  4. b =
  5.     0.0201    0.0402    0.0201
  6. a =
  7.     1.0000   -1.5610    0.6414

Further experiments

I tried to clone ruohoruotsi/Butterworth-Filter-Design.git and test with:

  1. #include "Butterworth.h"
  3. using namespace std;
  5. int main() {
  6.   float fs = 20000;
  7.   float fc = 500;
  8.   int order = 2;
  10.   vector<Biquad> coeffs; // second-order sections (sos)
  11.   Butterworth butterworth;
  13.   bool designedCorrectly = butterworth.loPass(fs, fc, 0, order, coeffs, 1.0f);
  15.   std::cout << "Designed correctly? " << designedCorrectly << std::endl;
  16.   std::cout << "a = [" << 1.f << ", " << coeffs[0].a1 << ", " << coeffs[0].a2
  17.             << "]" << std::endl;
  18.   std::cout << "b = [" << coeffs[0].b0 << ", " << coeffs[0].b1 << ", "
  19.             << coeffs[0].b2 << "]" << std::endl;
  20. }

It gives me:

  1. a = [1, 1.77863, -0.800803]
  2. b = [1, 2, 1]

Which again doesn't really fit.

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

答案1

得分: 1

以下是代码的一部分：

如果在z变换之后的传递函数定义为

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

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

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

  1. #include <iostream>
  2. #include <iomanip>
  3. #include <complex>
  4. #include <vector>
  5. #include <cmath>
  6. using namespace std;
  8. const double PI = 4.0 * atan( 1.0 );
  9. const complex<double> iPI( 0.0, PI );
  10. const double SQRT2 = sqrt( 2.0 );
  12. // const double factor = 0.5;    // use for Matlab comparison
  13. const double factor = 1.0;       // as original post
  15. //----------------------------------------------------------------------
  17. void getCoefficients( double fd, double fs, vector<double> &a, vector<double> &b )
  18. {
  19.    double alpha = 1.0 / tan( factor * PI * fd / fs );
  20.    b[0] = 1.0 / ( 1.0 + SQRT2 * alpha + alpha * alpha );
  21.    b[1] = 2 * b[0];
  22.    b[2] = b[0];
  23.    a[0] = 1.0;
  24.    a[1] = 2.0 * ( 1.0 - alpha * alpha ) * b[0];
  25.    a[2] = ( 1.0 - SQRT2 * alpha + alpha * alpha ) * b[0];
  26. }
  28. //----------------------------------------------------------------------
  30. void bode( const vector<double> &a, const vector<double> &b, double f, double fs, double &mag, double &phase )
  31. {
  32.    complex<double> invz = ( 1.0 - iPI * factor * f / fs ) / ( 1.0 + iPI * factor * f / fs );
  33.    complex<double> H = ( b[0] + b[1] * invz + b[2] * invz * invz ) / ( a[0] + a[1] * invz + a[2] * invz * invz );
  34.    mag = abs( H );
  35.    phase = atan2( H.imag(), H.real() );
  36. }
  38. //----------------------------------------------------------------------
  40. int main()
  41. {
  42.    vector<double> a(3), b(3);
  43.    double fs = 10000;
  44.    double fd = 1000;
  46.    // Get Butterworth coefficients in H(z) = ( b0 + b1/z + b2/z^2 ) / ( a0 + a1/z + a2/z^2 )
  47.    getCoefficients( fd, fs, a, b );
  48.    cout << "a:\t";
  49.    for ( double e : a ) cout << e << '\t';
  50.    cout << "\nb:\t";
  51.    for ( double e : b ) cout << e << '\t';
  52.    cout << "\n\n";
  54.    double logmin = 1.0, logmax = 4.0;
  55.    int nw = 31;
  56.    double dlog = ( logmax - logmin ) / ( nw - 1 );
  57.    #define FMT << scientific << setw(15) << setprecision(4) <<
  59.    cout FMT "f(Hz)" FMT "mag" FMT "db" FMT "phase(deg)" << '\n';
  60.    for ( int i = 0; i < nw; i++ )
  61.    {
  62.       double f, mag, phase;
  63.       f = pow( 10.0, logmin + i * dlog );
  64.       bode( a, b, f, fs, mag, phase );
  66.       double phasedeg = phase * 180.0 / PI;
  67.       double db = 20.0 * log( mag ) / log( 10.0 );
  68.       cout FMT f FMT mag FMT db FMT phasedeg << '\n';
  69.    }
  70. }

希望这有所帮助。如果你需要更多翻译，请提供具体的部分。

英文:

If the transfer function after a z-transform is defined as

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

then I think that the formulation for the coefficients a1 and a2 is producing values of the wrong sign. (At least, that's what I found when reworking it, but I might be wrong).

The following code gives the opposite signs for a1 and a2.

  1. #include &lt;iostream&gt;
  2. #include &lt;iomanip&gt;
  3. #include &lt;complex&gt;
  4. #include &lt;vector&gt;
  5. #include &lt;cmath&gt;
  6. using namespace std;
  7. const double PI = 4.0 * atan( 1.0 );
  8. const complex&lt;double&gt; iPI( 0.0, PI );
  9. const double SQRT2 = sqrt( 2.0 );
  10. // const double factor = 0.5;    // use for Matlab comparison
  11. const double factor = 1.0;       // as original post
  12. //----------------------------------------------------------------------
  13. void getCoefficients( double fd, double fs, vector&lt;double&gt; &amp;a, vector&lt;double&gt; &amp;b )
  14. {
  15. double alpha = 1.0 / tan( factor * PI * fd / fs );
  16. b[0] = 1.0 / ( 1.0 + SQRT2 * alpha + alpha * alpha );
  17. b[1] = 2 * b[0];
  18. b[2] = b[0];
  19. a[0] = 1.0;
  20. a[1] = 2.0 * ( 1.0 - alpha * alpha ) * b[0];
  21. a[2] = ( 1.0 - SQRT2 * alpha + alpha * alpha ) * b[0];
  22. }
  23. //----------------------------------------------------------------------
  24. void bode( const vector&lt;double&gt; &amp;a, const vector&lt;double&gt; &amp;b, double f, double fs, double &amp;mag, double &amp;phase )
  25. {
  26. complex&lt;double&gt; invz = ( 1.0 - iPI * factor * f / fs ) / ( 1.0 + iPI * factor * f / fs );
  27. complex&lt;double&gt; H = ( b[0] + b[1] * invz + b[2] * invz * invz ) / ( a[0] + a[1] * invz + a[2] * invz * invz );
  28. mag = abs( H );
  29. phase = atan2( H.imag(), H.real() );
  30. }
  31. //----------------------------------------------------------------------
  32. int main()
  33. {
  34. vector&lt;double&gt; a(3), b(3);
  35. double fs = 10000;
  36. double fd = 1000;
  37. // Get Butterworth coefficients in H(z) = ( b0 + b1/z + b2/z^2 ) / ( a0 + a1/z + a2/z^2 )
  38. getCoefficients( fd, fs, a, b );
  39. cout &lt;&lt; &quot;a:\t&quot;;
  40. for ( double e : a ) cout &lt;&lt; e &lt;&lt; &#39;\t&#39;;
  41. cout &lt;&lt; &quot;\nb:\t&quot;;
  42. for ( double e : b ) cout &lt;&lt; e &lt;&lt; &#39;\t&#39;;
  43. cout &lt;&lt; &quot;\n\n&quot;;
  44. double logmin = 1.0, logmax = 4.0;
  45. int nw = 31;
  46. double dlog = ( logmax - logmin ) / ( nw - 1 );
  47. #define FMT &lt;&lt; scientific &lt;&lt; setw(15) &lt;&lt; setprecision(4) &lt;&lt;
  48. cout FMT &quot;f(Hz)&quot; FMT &quot;mag&quot; FMT &quot;db&quot; FMT &quot;phase(deg)&quot; &lt;&lt; &#39;\n&#39;;
  49. for ( int i = 0; i &lt; nw; i++ )
  50. {
  51. double f, mag, phase;
  52. f = pow( 10.0, logmin + i * dlog );
  53. bode( a, b, f, fs, mag, phase );
  54. double phasedeg = phase * 180.0 / PI;
  55. double db = 20.0 * log( mag ) / log( 10.0 );
  56. cout FMT f FMT mag FMT db FMT phasedeg &lt;&lt; &#39;\n&#39;;
  57. }
  58. }

Output:

  1. a:	1	-1.14298	0.412802	
  2. b:	0.0674553	0.134911	0.0674553	
  3. f(Hz)            mag             db     phase(deg)
  4. 1.0000e+01     1.0000e+00    -3.7956e-08    -7.8347e-01
  5. 1.2589e+01     1.0000e+00    -9.5341e-08    -9.8635e-01
  6. 1.5849e+01     1.0000e+00    -2.3949e-07    -1.2418e+00
  7. 1.9953e+01     1.0000e+00    -6.0156e-07    -1.5634e+00
  8. 2.5119e+01     1.0000e+00    -1.5111e-06    -1.9683e+00
  9. 3.1623e+01     1.0000e+00    -3.7956e-06    -2.4783e+00
  10. 3.9811e+01     1.0000e+00    -9.5341e-06    -3.1205e+00
  11. 5.0119e+01     1.0000e+00    -2.3949e-05    -3.9296e+00
  12. 6.3096e+01     9.9999e-01    -6.0156e-05    -4.9494e+00
  13. 7.9433e+01     9.9998e-01    -1.5110e-04    -6.2354e+00
  14. 1.0000e+02     9.9996e-01    -3.7954e-04    -7.8588e+00
  15. 1.2589e+02     9.9989e-01    -9.5331e-04    -9.9113e+00
  16. 1.5849e+02     9.9972e-01    -2.3942e-03    -1.2513e+01
  17. 1.9953e+02     9.9931e-01    -6.0114e-03    -1.5821e+01
  18. 2.5119e+02     9.9826e-01    -1.5084e-02    -2.0052e+01
  19. 3.1623e+02     9.9566e-01    -3.7791e-02    -2.5501e+01
  20. 3.9811e+02     9.8920e-01    -9.4310e-02    -3.2581e+01
  21. 5.0119e+02     9.7352e-01    -2.3312e-01    -4.1849e+01
  22. 6.3096e+02     9.3720e-01    -5.6339e-01    -5.3957e+01
  23. 7.9433e+02     8.6132e-01    -1.2967e+00    -6.9313e+01
  24. 1.0000e+03     7.3050e-01    -2.7276e+00    -8.7273e+01
  25. 1.2589e+03     5.5943e-01    -5.0451e+00    -1.0563e+02
  26. 1.5849e+03     3.9180e-01    -8.1387e+00    -1.2189e+02
  27. 1.9953e+03     2.5949e-01    -1.1718e+01    -1.3493e+02
  28. 2.5119e+03     1.6715e-01    -1.5538e+01    -1.4496e+02
  29. 3.1623e+03     1.0636e-01    -1.9464e+01    -1.5262e+02
  30. 3.9811e+03     6.7339e-02    -2.3435e+01    -1.5850e+02
  31. 5.0119e+03     4.2546e-02    -2.7423e+01    -1.6305e+02
  32. 6.3096e+03     2.6859e-02    -3.1418e+01    -1.6660e+02
  33. 7.9433e+03     1.6951e-02    -3.5416e+01    -1.6939e+02
  34. 1.0000e+04     1.0696e-02    -3.9415e+01    -1.7159e+02

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

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

As for Matlab ... well, if you change the factor at the top of that code to the commented-out value (0.5) then you would get the Matlab values:

  1. a:      1       -1.56102        0.641352        
  2. b:      0.0200834       0.0401667       0.0200834

答案2

得分: 0

butter(2, 0.2) = [0.06745527, 0.13491055, 0.06745527], [1., -1.1429805, 0.4128016]

butter(2, 0.1) = [0.02008337, 0.04016673, 0.02008337], [1., -1.56101808, 0.64135154]

FF 值相对于 Matlab 的值差了一个 "2" 的因素。Butterworth.h 代码通过乘以 (1.0 + q*ita + ita*ita) 对 "b" 值进行了归一化。 "a" 值相对于 Matlab、Scipy 和原始 SO 帖子使用了负号，我没有检查它对结果的影响。

英文:

I dont have access to Matlab, but on python

butter(2,.2) = [0.06745527, 0.13491055, 0.06745527], [ 1., -1.1429805, 0.4128016]

and

butter(2,.1) = [0.02008337, 0.04016673, 0.02008337],[ 1. , -1.56101808, 0.64135154]

So your values from Matlab are off by a factor of "2" on FF (I once made same mistake and posted a comment in https://stackoverflow.com/questions/20924868)

The "Butterworth.h" code has normalized the "b" values by multiplying by (1.0 + q*ita + ita*ita). Your "a" values have minus signs compared to what Matlab, Scipy, and the original SO post use, I didnt check how that contributes too.

答案3

得分: 0

以下是翻译好的内容：

对于 Butterworth 滤波，您可以使用以下代码：

  1. void butter(float fc, float fs, bool lowPass = true) {
  2.     const float k = tanf(M_PI * fc / (fs * 2));
  3.     const float k2 = k * k;
  4.     const float q = 1 / sqrtf(2);
  5.     const float norm = 1 / (1 + k / q + k2);
  6.     if (lowPass) {
  7.         b0 = k2 * norm;
  8.         b1 = 2 * b0;
  9.         b2 = b0;
  10.     } else {
  11.         b0 = 1 * norm;
  12.         b1 = -2 * b0;
  13.         b2 = b0;
  14.     }
  15.     a1 = 2 * (k2 - 1) * norm;
  16.     a2 = (1 - k / q + k2) * norm;
  17. }

在 Matlab 中，您可以使用 butter(2, fc/fs) 获得完全相同的系数。

英文:

For Butterworth you could use this:

  1. void butter(float fc, float fs, bool lowPass = true) {
  2. const float k = tanf(M_PI * fc/(fs*2));
  3. const float k2 = k * k;
  4. const float q = 1 / sqrtf(2);
  5. const float norm = 1 / (1 + k / q + k2);
  6. if (lowPass) {
  7. b0 = k2 * norm;
  8. b1 = 2 * b0;
  9. b2 = b0;
  10. } else {
  11. b0 = 1 * norm;
  12. b1 = -2 * b0;
  13. b2 = b0;
  14. }
  15. a1 = 2 * (k2 - 1) * norm;
  16. a2 = (1 - k / q + k2) * norm;
  17. }

In Matlab you get the exact same coefficients with butter(2, fc/fs)

huangapple
  • 本文由 发表于 2023年5月7日 18:04:40
  • 转载请务必保留本文链接：https://go.coder-hub.com/76193236.html
  • c++
  • digital-filter
  • matlab
  • signal-processing
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定