Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
C++ source code to compute discrete Fourier transform
/*********************************************
* 離散フーリエ変換 *
* f(t) = 2 * sin(4 * t) + 3 * cos(2 * t) *
* ( 0 <= t < 2 * pi ) *
*********************************************/
#include <iostream> // for cout
#include <math.h> // for sin(), cos()
#include <stdio.h> // for printf()
#define N 100 // 分割数
#define CSV_DFT "DFT.csv" // 出力ファイル (DFT)
#define CSV_IDFT "IDFT.csv" // 出力ファイル (IDFT)
using namespace std;
/*
* 計算クラス
*/
class Calc
{
double SRC_re[N]; //元データの実部
double SRC_im[N]; //元データの虚部
double DFT_re[N]; //DFTの実部
double DFT_im[N]; //DFTの虚部
double IDFT_re[N]; //IDFTの実部
double IDFT_im[N]; //IDFTの虚部
public:
void makeSourceData(); // 元データ作成
void executeDFT(); // 離散フーリエ変換
void executeIDFT(); // 逆離散フーリエ変換
private:
double calcTerm(int n, double x); //各項計算
};
/*
* 元データ作成
*/
void Calc::makeSourceData()
{
int i;
for (i = 0; i < N; i++) {
SRC_re[i] = 2 * sin(4 * (2 * M_PI / N) * i)
+ 3 * cos(2 * (2 * M_PI / N) * i);
SRC_im[i] = 0.0;
}
}
/*
* 離散フーリエ変換
*/
void Calc::executeDFT()
{
int k, n; // LOOPインデックス
FILE *pf; // ファイルポインタ
// 出力ファイルOPEN
pf = fopen(CSV_DFT, "w");
// ヘッダ出力 ( k, 角周波数, 元データ(実部), 元データ(虚部), DFT(実部), DFT(虚部) )
fprintf(pf, "k,f,x_re,x_im,X_re,X_im\n");
// 計算・結果出力
for (k = 0; k < N; k++) {
DFT_re[k] = 0.0;
DFT_im[k] = 0.0;
for (n = 0; n < N; n++) {
DFT_re[k] += SRC_re[n] * ( cos((2 * M_PI / N) * k * n))
+ SRC_im[n] * ( sin((2 * M_PI / N) * k * n));
DFT_im[k] += SRC_re[n] * (-sin((2 * M_PI / N) * k * n))
+ SRC_im[n] * ( cos((2 * M_PI / N) * k * n));
}
fprintf(pf, "%d,%lf,%lf,%lf,%lf,%lf\n",
k, (2 * M_PI / N) * k, SRC_re[k], SRC_im[k], DFT_re[k], DFT_im[k]);
}
// 出力ファイルCLOSE
fclose(pf);
}
/*
* 逆離散フーリエ変換
*/
void Calc::executeIDFT()
{
int k, n; // LOOPインデックス
FILE *pf; // ファイルポインタ
// 出力ファイルOPEN
pf = fopen(CSV_IDFT, "w");
// ヘッダ出力 ( k, 角周波数, DFT(実部), DFT(虚部), IDFT(実部), IDFT(虚部) )
fprintf(pf, "k,f,X_re,X_im,x_re,x_im\n");
// 計算・結果出力
for (n = 0; n < N; n++) {
IDFT_re[n] = 0.0;
IDFT_im[n] = 0.0;
for (k = 0; k < N; k++) {
IDFT_re[n] += DFT_re[k] * (cos((2 * M_PI / N) * k * n))
- DFT_im[k] * (sin((2 * M_PI / N) * k * n));
IDFT_im[n] += DFT_re[k] * (sin((2 * M_PI / N) * k * n))
+ DFT_im[k] * (cos((2 * M_PI / N) * k * n));
}
IDFT_re[n] /= N;
IDFT_im[n] /= N;
fprintf(pf, "%d,%lf,%lf,%lf,%lf,%lf\n",
n, (2 * M_PI / N) * n, DFT_re[n], DFT_im[n], IDFT_re[n], IDFT_im[n]);
}
// 出力ファイルCLOSE
fclose(pf);
}
/*
* メイン処理
*/
int main()
{
try
{
// 計算クラスインスタンス化
Calc objCalc;
// 元データ作成
objCalc.makeSourceData();
// 離散フーリエ変換
objCalc.executeDFT();
// 逆離散フーリエ変換
objCalc.executeIDFT();
}
catch (...) {
cout << "例外発生!" << endl;
return -1;
}
// 正常終了
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.