音声IIRフィルタ設計（C言語）

開発メモ: バターワース IIRデジタルフィルタ、係数算出関数

バターワースデジタル IIRフィルタ、係数算出のメモ
伝達関数H(s) →(双一次変換)→H(z) 係数展開

IIRフィルタ係数格納構造体

IIR.h

typedef struct
{
    double b0, a0;
    double b1, a1;
    double b2, a2;
    double b3, a3;
    double z[6];
} Coef3OrderType; // 3次IIR

typedef struct
{
    double b0, a0;
    double b1, a1;
    double b2, a2;
    double z[4];
} Coef2OrderType; // 2次IIR

typedef struct
{
    double b0, a0;
    double b1, a1;
    double z[2];
} Coef1OrderType; // 1次IIR

バターワース1次　LPF/HPF

IIR_1dim.c

// Create 1-order filter coef
// coef:係数構造体　fs:サンプリング周波数　fc:カットオフ周波数
// type: 0=LPF, 1=HPF
void CreateCoefOrder_1(Coef1OrderType *coef, double fs, double fc, int type)
{
    double PI = 3.14159265358979323846;
    double fa = 1.0 / (2.0*PI) * tan(PI*fc / fs);
    double pfc = 2.0*PI*fa;

    memset(coef, 0, sizeof(Coef1OrderType));

    if (type == 0) // LPF
    {
        coef->b0 = 1.0;
        coef->b1 = (pfc - 1.0) / (pfc + 1.0);
        coef->a0 = pfc / (pfc + 1.0);
        coef->a1 = pfc / (pfc + 1.0);
    }
    else // if(type == 1) // HPF
    {
        coef->b0 = 1.0;
        coef->b1 = (pfc - 1.0) / (pfc + 1.0);
        coef->a0 = 1.0 / (pfc + 1.0);
        coef->a1 = -1.0 / (pfc + 1.0);
    }
}

バターワース2次　LPF/HPF、 2次APF

IIR_2dim.c

// Create 2-order filter coef
// coef:係数構造体　fs:サンプリング周波数　fc:カットオフ周波数
// type: 0=LPF, 1=HPF, 2=APF
void CreateCoefOrder_2(Coef2OrderType *coef, double fs, double fc, int type)
{
    double PI = 3.14159265358979323846;
    double fa = 1.0 / (2.0*PI) * tan(PI*fc / fs);
    double pfc = 2.0*PI*fa;
    double RT2 = sqrt(2.0);

    memset(coef, 0, sizeof(Coef2OrderType));

    if (type == 0) // LPF
    {
        coef->b0 = 1.0;
        coef->b1 = (-2.0 + 2.0*pfc*pfc) / (1 + RT2*pfc + pfc*pfc);
        coef->b2 = (1.0 - RT2*pfc + pfc*pfc) / (1 + RT2*pfc + pfc*pfc);

        coef->a0 = pfc*pfc / (1 + RT2*pfc + pfc*pfc);
        coef->a1 = 2.0*pfc*pfc / (1 + RT2*pfc + pfc*pfc);
        coef->a2 = pfc*pfc / (1 + RT2*pfc + pfc*pfc);
    }
    else if(type == 1) // HPF
    {
        coef->b0 = 1.0;
        coef->b1 = (2.0*pfc*pfc - 2.0) / (pfc*pfc + RT2*pfc + 1.0);
        coef->b2 = (pfc*pfc - RT2*pfc + 1.0) / (pfc*pfc + RT2*pfc + 1.0);

        coef->a0 = 1.0 / (pfc*pfc + RT2*pfc + 1.0);
        coef->a1 = -2.0 / (pfc*pfc + RT2*pfc + 1.0);
        coef->a2 = 1.0 / (pfc*pfc + RT2*pfc + 1.0);
    }
    else // if(type == 1) // APF
    {
        double w0 = 2.0 * PI*fc / fs;
        double alpha = sin(w0) / 2.0;

        coef->b0 = 1.0;
        coef->b1 = -2 * cos(w0) / (1 + alpha);
        coef->b2 = (1 - alpha) / (1 + alpha);

        coef->a0 = (1 - alpha) / (1 + alpha);
        coef->a1 = -2 * cos(w0) / (1 + alpha);
        coef->a2 = (1 + alpha) / (1 + alpha);
    }
}

バターワース　3次　LPF/HPF

IIR_2dim.c

// Create 3-order filter coef
// coef:係数構造体　fs:サンプリング周波数　fc:カットオフ周波数
// type: 0=LPF, 1=HPF
void CreateCoefOrder_3(Coef3OrderType *coef, double fs, double fc, int type)
{
    double PI = 3.14159265358979323846;
    double fa = 1.0 / (2.0*PI) * tan(PI*fc / fs);
    double pfc = 2.0*PI*fa;

    memset(coef, 0, sizeof(Coef3OrderType));

    if (type == 0) // LPF
    {
        double tmp = (1.0 + 2.0*pfc + 2.0*pfc*pfc + pfc*pfc*pfc);

        coef->b0 = 1.0;
        coef->b1 = (-3.0 - 2.0*pfc + 2.0*pfc*pfc + 3.0*pfc*pfc*pfc) / tmp;
        coef->b2 = (3.0 - 2.0*pfc - 2.0*pfc*pfc + 3.0*pfc*pfc*pfc) / tmp;
        coef->b3 = (-1.0 + 2.0*pfc - 2.0*pfc*pfc + pfc*pfc*pfc) / tmp;

        coef->a0 = (pfc*pfc*pfc) / tmp;
        coef->a1 = 3.0*(pfc*pfc*pfc) / tmp;
        coef->a2 = 3.0*(pfc*pfc*pfc) / tmp;
        coef->a3 = (pfc*pfc*pfc) / tmp;
    }
    else // if(type == 1) // HPF
    {
        double tmp = (1.0 + 2.0*pfc + 2.0*pfc*pfc + pfc*pfc*pfc);

        coef->b0 = 1.0;
        coef->b1 = (-3.0 - 2.0*pfc + 2.0*pfc*pfc + 3.0*pfc*pfc*pfc) / tmp;
        coef->b2 = (3.0 - 2.0*pfc - 2.0*pfc*pfc + 3.0*pfc*pfc*pfc) / tmp;
        coef->b3 = (-1.0 + 2.0*pfc - 2.0*pfc*pfc + pfc*pfc*pfc) / tmp;

        coef->a0 = 1.0 / tmp;
        coef->a1 = -3.0 / tmp;
        coef->a2 = 3.0 / tmp;
        coef->a3 = -1.0 / tmp;
    }
}

2次 ピーキングフィルタ（パラメトリックEQ用）

IIR_3dim.c

// coef:係数構造体　fs:サンプリング周波数　fc:カットオフ周波数 dBgain:　ゲイン(db) Q: 鋭度
void CreateCoefOrder2_PEAKING(Coef3OrderType *coef, double fs, double fc, double dBgain, double Q)
{
    double norm;
    double A, omega, sinval, cosval, alpha;
    // intermediate
    A = pow(10.0, (dBgain / 40));
    omega = 2.0*PI*fc / fs;
    sinval = sin(omega);
    cosval = cos(omega);
    alpha = sinval / (2.0*Q);

    memset(coef, 0, sizeof(Coef2OrderType));

    // result peaking
    coef->a0 = 1 + alpha *A;
    coef->a1 = -2.0 * cosval;
    coef->a2 = 1 - alpha * A;
    coef->b0 = 1 + alpha / A;
    coef->b1 = -2.0 * cosval;
    coef->b2 = 1 - alpha / A;

    norm = coef->b0;

    coef->b0 /= norm;
    coef->b1 /= norm;
    coef->b2 /= norm;
}

IIRフィルタ実行関数（1次、２次、３次）

引数：　cf:係数構造体　inp:入力
戻り値：　出力

IIR_dofil.c

//　1次IIR実行関数
double DoFilter1( Coef1OrderType *cf, double inp )
{
    double outp;

    // calc
    outp = inp*cf->a0 + cf->z[0]*cf->a1  - cf->z[1]*cf->b1;
    // unit delay
    cf->z[0] = inp;
    cf->z[1] = outp;

    return outp;
}

//　2次IIR実行関数
double DoFilter2( Coef2OrderType *cf, double inp )
{
    double outp;

    // calc
    outp = inp*cf->a0 + cf->z[0]*cf->a1 + cf->z[1]*cf->a2 
            - cf->z[2]*cf->b1 - cf->z[3]*cf->b2;
    // unit delay
    cf->z[1] = cf->z[0];
    cf->z[0] = inp;
    cf->z[3] = cf->z[2];
    cf->z[2] = outp;

    return outp;
}

//　3次IIR実行関数
double DoFilter3( Coef3OrderType *cf, double inp )
{
    double outp;

    // calc
    outp = inp*cf->a0 + cf->z[0]*cf->a1 + cf->z[1]*cf->a2 + cf->z[2]*cf->a3 
            - cf->z[3]*cf->b1 - cf->z[4]*cf->b2 - cf->z[5]*cf->b3;
    // unit delay
    cf->z[2] = cf->z[1];
    cf->z[1] = cf->z[0];
    cf->z[0] = inp;
    cf->z[5] = cf->z[4];
    cf->z[4] = cf->z[3];
    cf->z[3] = outp;

    return outp;
}