離散フーリエ変換

(DFT:Discrete Fourier Transform)

離散フーリエ変換 (DFT) は下記のように定義されています。

離散フーリエ変換(DFT:Discrete Fourier Transform)

    　  　N-1　　   jk                -2πi/N    
    A  =  Σ   a   W   ，      W  =  e   　         
    　k 　j=0　 j   N  　   　  N 　 　   　         


離散フーリエ逆変換(IDFT:Inverse Discrete Fourier Transform)

    　  　 1  N-1        -jk    
    a  =  --- Σ   A   W       
     k 　  N  j=0   j   N      


オイラーの公式(Euler's Formula)
  iθ
 e   ＝ cosθ＋i sinθ
    

オイラーの公式(Euler's Formula)を利用して、離散フーリエ変換、及び逆変換を簡単に実現できます。

//離散フーリエ変換(DFT:Discrete Fourier Transform)
void DFT(double ar[], double ai[], double AR[], double AI[], int Num)
{
	double theta = - 8.0 * atan(1.0) / Num;
	for (int k = 0; k < Num; k ++)
	{
		AR[k] = 0.0;
		AI[k] = 0.0;
		for (int j = 0; j < Num; j ++)
		{
			AR[k] += ar[j] * cos(j * k * theta) - ai[j] * sin(j * k * theta);
			AI[k] += ar[j] * sin(j * k * theta) + ai[j] * cos(j * k * theta);
		}
	}
}

//離散フーリエ逆変換(IDFT:Inverse Discrete Fourier Transform)
void IDFT(double AR[], double AI[], double ar[], double ai[], int Num)
{
	double theta = 8.0 * atan(1.0) / Num;
	for (int k = 0; k < Num; k ++)
	{
		ar[k] = 0.0;
		ai[k] = 0.0;
		for (int j = 0; j < Num; j ++)
		{
			ar[k] += AR[j] * cos(j * k * theta) - AI[j] * sin(j * k * theta);
			ai[k] += AR[j] * sin(j * k * theta) + AI[j] * cos(j * k * theta);
		}
		ar[k] /= Num;
		ai[k] /= Num;
	}
}

テスト：波形を配列に入れて、離散フーリエ変換と逆変換して、結果を画面に表示します。

////////////////////////////////////////////////////////
// テスト
//データを1000個にする
#define N 1000		//サンプリング周波数の1000Hzにする
double ar[N];		//データの実数部
double ai[N];		//データの虚数部
double AR[N];		//DFT後の実数部
double AI[N];		//DFT後の虚数部
double ar2[N];		//データの実数部
double ai2[N];		//データの虚数部

//初期データを設定する
void InitData()
{
	double alfa = 8.0 * atan(1.0) / N;
	for (int i = 0; i < N; i ++)
	{
		ar[i] = 10.0 * cos(i * alfa * 10)		//10HZの余弦波
				+10.0 * cos(i * alfa * 20)		//20HZの余弦波
				+10.0 * sin(i * alfa * 30);		//30HZの正弦波
		ai[i] = 0.0;
	}
}

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

	InitData();
	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(i, 100 - (int)ar[i]);
		else
			pDC->LineTo(i, 100 - (int)ar[i]);
	}
	DFT(ar, ai, AR, AI, N);

	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(i, 200 - (int)(AR[i] * 2 / N));
		else
			pDC->LineTo(i, 200 - (int)(AR[i] * 2 / N));
	}
	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(i, 300 - (int)(AI[i] * 2 / N));
		else
			pDC->LineTo(i, 300 - (int)(AI[i] * 2 / N));
	}

	IDFT(AR, AI, ar2, ai2, N);

	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(i, 400 - (int)ar2[i]);
		else
			pDC->LineTo(i, 400 - (int)ar2[i]);
	}
	for (int i = 0; i < N; i ++)
	{
		if (i == 0)
			pDC->MoveTo(i, 500 - (int)ai2[i]);
		else
			pDC->LineTo(i, 500 - (int)ai2[i]);
	}
}