c - 如何计算双二阶滤波器系数？

Question

这是我用于计算双二阶滤波器系数的 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;
}

Answer 1

代码未完全实现引用的Cookbook formulae for audio equalizer biquad filter coefficients等式#4。

有 3 种情况，仅代码代码情况 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 数学进行了所有计算，建议重新编码以仅使用双。 float *coeff --> double *coeff。