離散フーリエ逆変換を利用で、LPF-FIR係数作成 LPFフィルタ係数の周波数特性を持つ配列を離散フーリエ逆変換 (IDFT)して、FIRの係数を作る手法のメモです。 １．LPFフィルタ係数の周波数特性を周波数特性を持つ配列にセットします。虚数部をゼロにして、実数部は、カットオフ周波数未満のデータを１にして、カットオフ周波数以上のデータをゼロにします。当然、 周波数特性は左右対称ですので、右半分を忘れずに。 例：サンプル周波数＝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.00000000LPFのテスト。周波数の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);
}
|
