素朴なアルゴリズムで離散フーリエ変換(DFT)を計算するC++プログラム

離散フーリエ変換(DFT:Discrete Fourier Transform)とは、関数を等間隔でサンプリングして得られた有限個の標本列を、複素正弦波の有限な線形結合における係数列へと変換する手法です。求められた係数は周波数の順に並べられ、元の標本値と同じ情報を持つため、サンプリングされた関数をその元の領域(多くの場合、時間や直線上の位置)から周波数領域へと変換できます。

DFTの基本式

N点の標本列 x(0), x(1), …, x(N−1) に対して、k番目のDFT係数は次のように定義されます。

X(k) = Σ x(i) × e^−2πik/N(i = 0 ～ N−1)

オイラーの公式により、指数部は余弦(コサイン)と正弦(サイン)の項に分離できるため、実部と虚部はそれぞれ以下のように計算できます。

• 実部:Re X(k) = Σ x(i) × cos(2πik / N)
• 虚部:Im X(k) = Σ x(i) × sin(2πik / N)

アルゴリズム

開始
    変数 M を宣言し、適当な整数で初期化する
    配列 function[M] を宣言する
    i = 0 から M−1 まで繰り返す
        function[i] = (((a * (double)i) + (b * (double)i)) − c)
    繰り返し終了
    配列 sine[M]、cosine[M] を宣言する
    i = 0 から M−1 まで繰り返す
        cosine[i] = cos((2 * i * k * PI) / M)
        sine[i] = sin((2 * i * k * PI) / M)
    繰り返し終了
    DFT_Coeff 型の配列 dft_value[k] を宣言する
    j = 0 から k−1 まで繰り返す
        i = 0 から M−1 まで繰り返す
            dft_value[j].real += function[i] * cosine[i]
            dft_value[j].img += function[i] * sine[i]
        繰り返し終了
    繰り返し終了
    結果を出力する
終了

C++による実装例

#include<iostream>
#include<math.h>
using namespace std;
#define PI 3.14159265

class DFT_Coeff {
public:
    double real, img;
    DFT_Coeff() {
        real = 0.0;
        img = 0.0;
    }
};

int main(int argc, char **argv) {
    int M = 10;
    cout << "Enter the coefficient of simple linear function:\n";
    cout << "ax + by = c\n";
    double a, b, c;
    cin >> a >> b >> c;

    double function[M];
    for (int i = 0; i < M; i++) {
        function[i] = (((a * (double)i) + (b * (double)i)) - c);
    }

    cout << "Enter the max K value: ";
    int k;
    cin >> k;

    double cosine[M];
    double sine[M];
    for (int i = 0; i < M; i++) {
        cosine[i] = cos((2 * i * k * PI) / M);
        sine[i] = sin((2 * i * k * PI) / M);
    }

    DFT_Coeff dft_value[k];
    cout << "The coefficients are: \n";
    for (int j = 0; j < k; j++) {
        for (int i = 0; i < M; i++) {
            dft_value[j].real += function[i] * cosine[i];
            dft_value[j].img += function[i] * sine[i];
        }
        cout << "(" << dft_value[j].real << ") - " << "(" << dft_value[j].img << " i)\n";
    }
}

実行結果

Enter the coefficient of simple linear function:
ax + by = c
4 5 6
Enter the max K value: 10
The coefficients are:
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)
(345) - (-1.64772e-05 i)

補足:計算量と注意点

この素朴な(naive)アプローチでは、k個の係数それぞれについてM回の内積計算を行うため、時間計算量は O(M×k) となります。データサイズが大きくなる場合は、O(N log N) で計算できる FFT(高速フーリエ変換)の利用が推奨されます。一方、小規模なデータやDFTの仕組みを学習する目的には、このシンプルな実装が非常に分かりやすく適しています。

また、上記のコードでは三角関数の引数が固定の k に依存しているため、すべての j に対して同じ係数値が出力されます。厳密なDFTを実装する場合は、cos・sin の中の周波数インデックスを外側のループ変数 j に対応させるように修正してください。