C++/STLでDFT、IDFTしてみる。

C++/STLでモダンに複素DFT、IDFTを実装。なんか、もやもやしていたけどすっきりした。やってみたら簡単だった。


#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>

#include <iostream>
#include <complex>
#include <cmath>	// sqrt, sin, cos
#include <iomanip>	// setprecision

const double pi = 3.1415926535897931f;

const int length = 32;

// initialize array
template <class T> 
void inita(std::complex<T> *data, int length) {
	for(int i=0; i<length; i++) {
		data[i] = std::complex<T>(0.5f+std::sin(2*pi*i/length+4), std::cos(4*2*pi*i/length));
		//data[i] = std::complex<double>(2, 3);
	}
}

// display array
template <class T>
void dispa(std::complex<T> *data, int length) {
	for(int i=0; i<length; i++) {
		std::wcout << i << L":" << std::fixed << std::setprecision(5) << data[i] << L" " << abs(data[i]) << std::endl;
	}
}

// deep copy
template <class T>
void copya(std::complex<T> *data, int length, std::complex<T> *dest) {
	for(int i=0; i<length; i++) {
		dest[i] = data[i];
	}
}

// DFT
template <class T>
void worka(std::complex<T> *data, int length, std::complex<T> *res) {
	for(int j=0; j<length; j++) {
		res[j] = std::complex<T>(0.0f, 0.0f);
		for(int i=0; i<length; i++) {
			res[j] += exp(std::complex<T>(0,  2*pi*i*j/length)) * data[i];
		}
	}
}

// IDFT
template <class T>
void workb(std::complex<T> *data, int length, std::complex<T> *res) {
	for(int j=0; j<length; j++) {
		res[j] = std::complex<T>(0.0f, 0.0f);
		for(int i=0; i<length; i++) {
			res[j] += exp(std::complex<T>(0, -2*pi*i*j/length)) * data[i];
		}
		res[j] /= length;
	}
}

void main(void) {
	_CrtSetReportMode(_CRT_WARN,  _CRTDBG_MODE_FILE);
	_CrtSetReportMode(_CRT_ERROR, _CRTDBG_MODE_FILE | _CRTDBG_MODE_DEBUG);
	_CrtSetReportFile(_CRT_WARN,  _CRTDBG_FILE_STDERR);
	_CrtSetReportFile(_CRT_ERROR, _CRTDBG_FILE_STDERR);
	//#define d_new new(_NORMAL_BLOCK, __FILE__, __LINE__)

	using atype = double;
	std::complex<atype> *data, *res;
	data = new std::complex<atype>[length];
	res  = new std::complex<atype>[length];
	
	std::wcout << L"-- orig --" << std::endl;
	inita(data, length);
	dispa(data, length);
	
	std::wcout << L"-- DFT --" << std::endl;
	worka(data, length, res);
	copya(res, length, data);
	dispa(data, length);
	
	std::wcout << L"-- IDFT --" << std::endl;
	workb(data, length, res);
	dispa(res, length);

	delete [] data;
	delete [] res;
	
	_CrtDumpMemoryLeaks();

	return;
}

ビルドはVC++で cl /EHa complex.cpp

十数年前に学校で履修してわけがわからないまま単位取得したテイラー展開フーリエ展開ラプラス変換、それと線形代数の行列・逆行列はすべて関連しているんだと、なぜか未明に気づいてはっとした。数学が得意な人からしたら当たり前じゃんてなるだろうが、数学が嫌いな人からしたらそれぞれは離散的で結びつかない。こういうことはもっと早く教えてほしかった。ウィキペディアの説明ではとても理解できない。



忘れないうちにコードのメモ。

main()内の_CrtSetReportModeあたりはCランタイム組み込みのメモリリーク検出デバッグ機能の有効化。最後の _CrtDumpMemoryLeaks();を実行するときに開放されていないメモリがあれば警告してくる。#define d_newを有効にしてヒープのnewをd_newに変更すればどこでメモリリークしたか表示する。

using atype = のところで全体の精度を指定する。どうせfloatかdoubleしか指定しないが、一応テンプレート対応てことで。Old C Styleのtypedef ~ と同じ。length = の指定で長さは任意に変更できる。FFTではないのでn^2長さでなくても構わない。

inita()で初期データを作成する。サンプルはY=0.5+SIN(2πt+4) + COS(4*2πt) i のオフセットのある実数部正弦波+虚数部余弦波のツートーン。

worka()、workb()がDFT、IDFTのメイン処理。同じデータ長の実数部と虚数部で90度位相差がある二つの正弦波を0倍(DC)~N-1倍の振動数まで順に掛け合わせる。DFTとIDFTは掛け合わせる振動の向きが正負逆になる。オイラーの公式で実数と虚数に分離してSIN、COSの組み合わせで計算すると逆行列になる。

プログラミング書籍なら翔泳社の通販『SEshop』

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト /  変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト /  変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト /  変更 )

%s と連携中