Go to Qiita Advent Calendar 2024 Top
LoginSignup
19
13

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@Micochan

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

Posted at

#開発メモ: バターワース 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;
}
19
13
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
19
13

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?