現在位置: Toshi Software Gallery > 開発メモ > IIRフィルタ係数算出

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

FabulouMP3に使用しているバターワースデジタル IIRフィルタ、係数算出。
伝達関数H(s) →(双一次変換)→H(z) 係数展開　（２次LPF/HPF計算ノート）

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

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/HP

// 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

// 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

// 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用）

// 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:入力
戻り値：　出力

//　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;
}