プログラミングメモ →目次

離散フーリエ逆変換を利用で、LPF-FIR係数作成

LPFフィルタ係数の周波数特性を持つ配列を離散フーリエ逆変換 (IDFT)して、FIRの係数を作る手法のメモです。 

1.LPFフィルタ係数の周波数特性を周波数特性を持つ配列にセットします。虚数部をゼロにして、実数部は、カットオフ周波数未満のデータを1にして、カットオフ周波数以上のデータをゼロにします。当然、 周波数特性は左右対称ですので、右半分を忘れずに。
2.周波数特性を持つ配列を逆離散フーリエ変換処理して、振幅/時間領域のデータを作ります。 
3.振幅/時間領域のデータは、両端に集中しているので、振幅/時間領域のデータを持つ配列を半周シフトします。両端のデータを真ん中に移動して、FIR係数を作ります。
4.窓関数を掛けて、FIR係数の段数を少なくします。 
5.作り立ての段数の少ないFIR係数を、一度に離散フーリエ変換して、周波数特性を確認します。

離散フーリエ逆変換、及び処理関数は、離散フーリエ変換を参照。

例:サンプル周波数=1000HZ、カットオフ周波数=100HZ、201段のLPF-FIRを作ります。

#define N 1000			//サンプリング周波数の1000Hzにする
double ar[N];			//データの実数部
double ai[N];			//データの虚数部
double AR[N];			//DFT後の実数部
double AI[N];			//DFT後の虚数部
double AR2[N];			//DFT後の実数部
double AI2[N];			//DFT後の虚数部
double ar_shift[N];		//シフトしたデータの実数部
double ar_window[N];	//シフトしたデータの実数部
#define M 201			//フィルタの段数
#define FC 100			//カットオフ周波数


void MakeSimpleLpfFirCoefficient(int iFreqCut)
{
	ZeroMemory(AI, sizeof(AI));
	ZeroMemory(AR, sizeof(AR));
	ZeroMemory(ar_shift, sizeof(ar_shift));
	ZeroMemory(ar_window, sizeof(ar_window));


	//周波数特性をセットする
	for (int i = 0; i < iFreqCut; i ++)
	{
		AR[i] = AR[N-1-i] = 1.0;
	}

	//離散フーリエ逆変換して、振幅/時間領域のデータを算出する
	IDFT(AR, AI, ar, ai, N);

	//シフトする
	for (int i = 0; i < N; i ++)
	{
		ar_shift[i] = ar[(i+(N/2)) % N];
	}

	//フィルタ係数に窓関数を掛ける //ハン窓
	
	for (int i = 0; i < M; i ++)
	{
		ar_window[i] = ar_shift[N/2 - M/2 + i];
		ar_window[i] = ar_window[i] * (0.5 - 0.5 * cos(8.0 * atan(1.0) * i / (M-1)));
	}

	//離散フーリエ変換して、フィルタ係数の周波数特性を確認する
	ZeroMemory(ai, sizeof(ai));
	DFT(ar_window, ai, AR2, AI2, N);
}

結果を確認するために、View描画を実装しています。

void CMakeLPF_Fir_CoefView::OnDraw(CDC* pDC)
{
	CMakeLPF_Fir_CoefDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (!pDoc)
		return;

	MakeSimpleLpfFirCoefficient(FC);

	pDC->TextOutW(10,30, CString("セットした周波数特性"));

	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(10+i, 100 - (int)(AR[i]*20));
		else
			pDC->LineTo(10+i, 100 - (int)(AR[i]*20));
	}

	pDC->TextOutW(10,130, CString("逆DFTして得たデータ"));
	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(10+i, 200 - (int)(ar[i]*N/4));
		else
			pDC->LineTo(10+i, 200 - (int)(ar[i]*N/4));
	}

	pDC->TextOutW(10,230, CString("シフトしたデータ→フィルタ係数"));
	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(10+i, 300 - (int)(ar_shift[i]*N/4));
		else
			pDC->LineTo(10+i, 300 - (int)(ar_shift[i]*N/4));
	}

	pDC->TextOutW(10,330, CString("窓を掛けたフィルタ係数→最終フィルタ係数"));
	for (int i = 0; i < M; i ++)
	{
		if (i == 0)
			pDC->MoveTo(10+i, 400 - (int)(ar_window[i]*N/4));
		else
			pDC->LineTo(10+i, 400 - (int)(ar_window[i]*N/4));
	}


	pDC->TextOutW(10,430, CString("最終フィルタ係数の周波数特性(確認用)"));
	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(10+i, 500 - (int)((AR2[i]*AR2[i]+AI2[i]*AI2[i])*20));
		else
			pDC->LineTo(10+i, 500 - (int)((AR2[i]*AR2[i]+AI2[i]*AI2[i])*20));
	}
}
 

LPF-FIRの係数は、下記のように作りました。

Coef[0]=-0.00000000
Coef[1]=-0.00000045
Coef[2]=-0.00000295
Coef[3]=-0.00000671
Coef[4]=-0.00000745
Coef[5]=-0.00000000
Coef[6]=0.00001711
Coef[7]=0.00003808
Coef[8]=0.00005024
Coef[9]=0.00003970
Coef[10]=-0.00000000
Coef[11]=-0.00006051
Coef[12]=-0.00011770
Coef[13]=-0.00013951
Coef[14]=-0.00010100
Coef[15]=-0.00000000
Coef[16]=0.00013454
Coef[17]=0.00024818
Coef[18]=0.00028097
Coef[19]=0.00019538
Coef[20]=0.00000000
Coef[21]=-0.00024337
Coef[22]=-0.00043639
Coef[23]=-0.00048161
Coef[24]=-0.00032725
Coef[25]=-0.00000000
Coef[26]=0.00039158
Coef[27]=0.00068991
Coef[28]=0.00074920
Coef[29]=0.00050154
Coef[30]=-0.00000000
Coef[31]=-0.00058439
Coef[32]=-0.00101747
Coef[33]=-0.00109274
Coef[34]=-0.00072401
Coef[35]=-0.00000000
Coef[36]=0.00082800
Coef[37]=0.00142946
Coef[38]=0.00152309
Coef[39]=0.00100166
Coef[40]=-0.00000000
Coef[41]=-0.00113009
Coef[42]=-0.00193898
Coef[43]=-0.00205401
Coef[44]=-0.00134346
Coef[45]=-0.00000000
Coef[46]=0.00150070
Coef[47]=0.00256323
Coef[48]=0.00270378
Coef[49]=0.00176145
Coef[50]=0.00000000
Coef[51]=-0.00195355
Coef[52]=-0.00332601
Coef[53]=-0.00349797
Coef[54]=-0.00227261
Coef[55]=0.00000000
Coef[56]=0.00250841
Coef[57]=0.00426190
Coef[58]=0.00447405
Coef[59]=0.00290209
Coef[60]=-0.00000000
Coef[61]=-0.00319513
Coef[62]=-0.00542374
Coef[63]=-0.00568990
Coef[64]=-0.00368919
Coef[65]=0.00000000
Coef[66]=0.00406144
Coef[67]=0.00689690
Coef[68]=0.00724023
Coef[69]=0.00469900
Coef[70]=0.00000000
Coef[71]=-0.00518855
Coef[72]=-0.00882892
Coef[73]=-0.00929117
Coef[74]=-0.00604758
Coef[75]=-0.00000000
Coef[76]=0.00672652
Coef[77]=0.01149761
Coef[78]=0.01216216
Coef[79]=0.00796305
Coef[80]=-0.00000000
Coef[81]=-0.00898518
Coef[82]=-0.01549276
Coef[83]=-0.01655200
Coef[84]=-0.01096116
Coef[85]=-0.00000000
Coef[86]=0.01271997
Coef[87]=0.02231694
Coef[88]=0.02433024
Coef[89]=0.01649955
Coef[90]=-0.00000000
Coef[91]=-0.02037050
Coef[92]=-0.03723907
Coef[93]=-0.04271959
Coef[94]=-0.03090315
Coef[95]=-0.00000000
Coef[96]=0.04658760
Coef[97]=0.10068333
Coef[98]=0.15121401
Coef[99]=0.18705108
Coef[100]=0.20000000	<-中心
Coef[101]=0.18705108
Coef[102]=0.15121401
Coef[103]=0.10068333
Coef[104]=0.04658760
Coef[105]=0.00000000
Coef[106]=-0.03090315
Coef[107]=-0.04271959
Coef[108]=-0.03723907
Coef[109]=-0.02037050
Coef[110]=0.00000000
Coef[111]=0.01649955
Coef[112]=0.02433024
Coef[113]=0.02231694
Coef[114]=0.01271997
Coef[115]=0.00000000
Coef[116]=-0.01096116
Coef[117]=-0.01655200
Coef[118]=-0.01549276
Coef[119]=-0.00898518
Coef[120]=0.00000000
Coef[121]=0.00796305
Coef[122]=0.01216216
Coef[123]=0.01149761
Coef[124]=0.00672652
Coef[125]=-0.00000000
Coef[126]=-0.00604758
Coef[127]=-0.00929117
Coef[128]=-0.00882892
Coef[129]=-0.00518855
Coef[130]=0.00000000
Coef[131]=0.00469900
Coef[132]=0.00724023
Coef[133]=0.00689690
Coef[134]=0.00406144
Coef[135]=0.00000000
Coef[136]=-0.00368919
Coef[137]=-0.00568990
Coef[138]=-0.00542374
Coef[139]=-0.00319513
Coef[140]=-0.00000000
Coef[141]=0.00290209
Coef[142]=0.00447405
Coef[143]=0.00426190
Coef[144]=0.00250841
Coef[145]=-0.00000000
Coef[146]=-0.00227261
Coef[147]=-0.00349797
Coef[148]=-0.00332601
Coef[149]=-0.00195355
Coef[150]=0.00000000
Coef[151]=0.00176145
Coef[152]=0.00270378
Coef[153]=0.00256323
Coef[154]=0.00150070
Coef[155]=0.00000000
Coef[156]=-0.00134346
Coef[157]=-0.00205401
Coef[158]=-0.00193898
Coef[159]=-0.00113009
Coef[160]=0.00000000
Coef[161]=0.00100166
Coef[162]=0.00152309
Coef[163]=0.00142946
Coef[164]=0.00082800
Coef[165]=0.00000000
Coef[166]=-0.00072401
Coef[167]=-0.00109274
Coef[168]=-0.00101747
Coef[169]=-0.00058439
Coef[170]=0.00000000
Coef[171]=0.00050154
Coef[172]=0.00074920
Coef[173]=0.00068991
Coef[174]=0.00039158
Coef[175]=0.00000000
Coef[176]=-0.00032725
Coef[177]=-0.00048161
Coef[178]=-0.00043639
Coef[179]=-0.00024337
Coef[180]=-0.00000000
Coef[181]=0.00019538
Coef[182]=0.00028097
Coef[183]=0.00024818
Coef[184]=0.00013454
Coef[185]=0.00000000
Coef[186]=-0.00010100
Coef[187]=-0.00013951
Coef[188]=-0.00011770
Coef[189]=-0.00006051
Coef[190]=-0.00000000
Coef[191]=0.00003970
Coef[192]=0.00005024
Coef[193]=0.00003808
Coef[194]=0.00001711
Coef[195]=0.00000000
Coef[196]=-0.00000745
Coef[197]=-0.00000671
Coef[198]=-0.00000295
Coef[199]=-0.00000045
Coef[200]=-0.00000000
LPFのテスト。周波数の25HZと200HZの正弦波の混じり合っている波形に、上記のLPFを掛けて、結果を見ます。

LPFのテスト用のプログラム:

	pDC->TextOutW(10,540, CString("テスト用波形"));
	double wav[1000];
	for (int i = 0; i < 1000; i ++)
	{
		wav[i] = 30.0 * sin(i * 8.0 * atan(1.0) * FC / 4 / N) + 10.0 * sin(i * 8.0 * atan(1.0) * FC * 2 / N);

		if (i == 0)
			pDC->MoveTo(10+i, 600 - (int)wav[i]);
		else
			pDC->LineTo(10+i, 600 - (int)wav[i]);
	}

	pDC->TextOutW(10,650, CString("LPFを掛けた後の波形"));
	for (int i = 0; i < 1000-M; i ++)
	{
		double dat = 0;
		for (int j = 0; j < M; j ++)
		{
			dat += wav[i+j] * ar_window[j];
		}

		if (i == 0)
			pDC->MoveTo(10+i, 700 - (int)dat);
		else
			pDC->LineTo(10+i, 700 - (int)dat);
	}