c - 如何计算双二阶滤波器系数?
问题描述
这是我用于计算双二阶滤波器系数的 C 代码。它们中的大多数是正确的,但不适用于峰值滤波器、低搁架和高搁架。
请帮忙看看我的代码是否有问题。
例如设置FS=48000,FC=1000,增益=10dB,Q=0.707,系数如下
Low Pass : 0.0039160769 0.00783215370 0.0039160769 -1.8153178692 0.8309822083
High Pass : 0.9115750194 -1.82315003870 0.9115750194 -1.8153178692 0.8309822083
Band Pass Peak Q : -0.0597477853 0.00000000000 -1.6619073153 -1.8153178692 0.8309822083
Band Pass Constant Q : 0.0845088884 0.00000000000 -0.0845088884 -1.8153178692 0.8309822083
Notch : 0.9154911041 -1.81531786920 0.9154911041 -1.8153178692 0.8309822083
All Pass : 0.8309822083 -1.81531786920 1.0000000000 -1.8153178692 0.8309822083
***Peaking EQ : 1.1067041159 -1.88503801820 0.7945997715 -1.8850380182 0.9013039470***
***Low Shelving : 1.0553481579 -1.85152983670 0.8246886730 -1.8612760305 0.8702906966***
***High Shelving : 2.9964308739 -5.57718515400 2.6077659130 -1.7544257641 0.7814375758***
预期的
The value for peak =
a0 = 1.18273168232815 a1 = -1.8153179156742152 a2 = 0.648250540080863 b1 = -1.8153179156742152 b2 = 0.8309822224090129
Low shelving
a0 = 1.074230188304032 a1 = -1.7984055763195854 a2 = 0.7737109074277086 b1 = -1.8153410827045682 b2 = 0.8310055893467576
High shelving
a0 = 3.0368659226630146 a1 = -5.757548057807517 a2 = 2.7363466417866915 b1 = -1.8153410827045682 b2 = 0.8310055893467576
比较的答案来自http://www.earlevel.com/main/2013/10/13/biquad-calculator-v2/
参考: http ://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
void audio_eq_coefficient(double Q, double dB, double FS, double FC,
int type, float *coeff) {
double n0, n1, n2, d1, d2;
double b0, b1, b2, a0, a1, a2;
double w0, cosW, sinW, A, alpha, beta;
// intermediate parameters
w0 = 2*M_PI*(FC/FS);
cosW = cos(w0);
sinW = sin(w0);
A = pow(10, (dB/40));
//printf("A = %lf, dB/40 = %lf\n", A, dB/40);
alpha = sinW/(2*Q);
beta = pow(A, 0.5)/Q;
switch (type) {
case 0: //Low Pass
b0 = (1 - cosW)/2;
b1 = 1 - cosW;
b2 = (1 - cosW)/2;
a0 = 1 + alpha;
a1 = -2*cosW;
a2 = 1 - alpha;
break;
case 1: //High Pass
b0 = (1 + cosW)/2;
b1 = -(1 + cosW);
b2 = (1 + cosW)/2;
a0 = (1 + alpha);
a1 = -2*cosW;
a2 = 1 - alpha;
break;
case 2: //Band Pass Peak gain Q
b0 = sinW/2;
b1 = 0;
b0 = -sinW/2;
a0 = 1 + alpha;
a1 = -2*cosW;
a2 = 1 - alpha;
break;
case 3: //Band Pass constant 0dB gain
b0 = alpha;
b1 = 0;
b2 = -alpha;
a0 = 1 + alpha;
a1 = -2*cosW;
a2 = 1 - alpha;
break;
case 4: //Notch Filter
b0 = 1;
b1 = -2*cosW;
b2 = 1;
a0 = 1 + alpha;
a1 = -2*cosW;
a2 = 1 - alpha;
break;
case 5: //All pass
b0 = 1 - alpha;
b1 = -2*cosW;
b2 = 1 + alpha;
a0 = 1 + alpha;
a1 = -2*cosW;
a2 = 1 - alpha;
break;
case 6: // Peaking EQ Filter
b0 = 1 + (alpha*A);
b1 = -2*cosW;
b2 = 1 - (alpha*A);
a0 = 1 + (alpha/A);
a1 = -2*cosW;
a2 = 1 - (alpha/A);
break;
case 7: //Low Shelving Filter
b0 = A*((A + 1) - ((A -1)*cosW) + beta*sinW);
b1 = 2*A*((A - 1) - ((A + 1)*cosW));
b2 = A*((A + 1) - (A - 1)*cosW - beta*sinW);
a0 = ((A + 1) + (A - 1)*cosW + beta*sinW);
a1 = -2*((A - 1) + (A + 1)*cosW);
a2 = ((A + 1) + (A - 1)*cosW - (beta*sinW));
break;
case 8: //High Shelving Filter
b0 = A*((A + 1) + (A - 1)*cosW + beta*sinW);
b1 = -2*A*((A - 1) + (A + 1)*cosW);
b2 = A*((A + 1) + (A - 1)*cosW -beta*sinW);
a0 = ((A + 1) - (A - 1)*cosW + beta*sinW);
a1 = 2*((A - 1) - (A + 1)*cosW);
a2 = ((A + 1) - (A - 1)*cosW - beta*sinW);
break;
default:
break;
};
n0 = b0/a0;
n1 = b1/a0;
n2 = b2/a0;
d1 = a1/a0;
d2 = a2/a0;
//printf("%lf %lf %lf %lf %lf\n", n0, n1, n2, d1, d2);
coeff[0] = (float)n0;
coeff[1] = (float)n1;
coeff[2] = (float)n2;
coeff[3] = (float)d1;
coeff[4] = (float)d2;
}
static char str_filter[9][30] = {"Low Pass", "High Pass", "Band Pass Peak Q", "Band Pass Constant Q",
"Notch", "All Pass", "Peaking EQ", "Low Shelving", "High Shelving"};
int main(){
float coeff[5];
int i;
for (i = 0; i < 9; i++) {
audio_eq_coefficient(0.707, 10, 48000, 1000, i, coeff);
printf("%-30s: %13.10f %13.10f0 %13.10f %13.10f %13.10f\n",
str_filter[i], coeff[0], coeff[1], coeff[2],
coeff[3], coeff[4]);
}
return 0;
}
解决方案
代码没有完全实现音频均衡器双二阶滤波器系数方程 #4的参考Cookbook 公式。
有 3 种情况,代码仅编码 case Q alpha = sinW/(2*Q);
。缺少的代码对应于过滤器 OP 有问题:峰值(案例:BW)和搁置(案例:S)。
case: Q
α=sin(ω0)/(2⋅Q)
case: BW
α=sin(ω0)⋅sinh(ln(2)/2⋅BW⋅ω0/sin(ω0))
case: S
α=sin(ω0)/22⋅√((A+1/A)⋅(1/S−1)+2)
注意:由于没有说明使用float
和编码这些方程的理由,至少已经用double
数学进行了所有计算,建议重新编码以仅使用double
. float *coeff
--> double *coeff
。
推荐阅读
- reactjs - 自定义钩子能否确定调用它的组件?
- angular - 如何使用角度从html中的数组中提取前两个字符?
- python - 如何每 x 秒重复调用 API 来更新坐标?
- javascript - 文件夹结构中的JS导入
- mysql - 如何在没有唯一值的两个表之间设置外键?
- reactjs - 仅当未使用 @azure/msal-react 进行身份验证时才重定向 onLoad
- javascript - 除非我刷新页面,否则 Vue data() 对象输入 v-model 不起作用
- postgresql - checkpoint_completion_target 被忽略
- python - 将两个不相等的列表压缩到字典中。如果 key 没有值,则不给出。忽略没有键的值
- vb.net - SQL Server 通配符 LIKE 搜索