英文:
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:
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);
}
};
Here the main:
int main()
{
const float fs = 10000;
Filter biquad(fs);
biquad.butterworthLowPass(1000);
biquad.displayCoefficients();
}
When taking the computed coefficients and using Matlab to display the filter response, I get a -15dB gain which is not what I expect.
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 gives me something very different:
>> fs = 10000;
>> fc = 1000;
>> [b, a] = butter(2, 1000/10000)
b =
0.0201 0.0402 0.0201
a =
1.0000 -1.5610 0.6414
Further experiments
I tried to clone ruohoruotsi/Butterworth-Filter-Design.git and test with:
#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;
}
It gives me:
a = [1, 1.77863, -0.800803]
b = [1, 2, 1]
Which again doesn't really fit.
答案1
得分: 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';
}
}
希望这有所帮助。如果你需要更多翻译，请提供具体的部分。
英文:
If the transfer function after a z-transform is defined as
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.
#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';
}
}
Output:
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
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:
a: 1 -1.56102 0.641352
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 滤波，您可以使用以下代码：
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 中，您可以使用 butter(2, fc/fs) 获得完全相同的系数。
英文:
For Butterworth you could use this:
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;
}
In Matlab you get the exact same coefficients with butter(2, fc/fs)
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论