離散フーリエ変換をC++で実装した

なぜかwolfram alphaの計算結果とずれてるけど、wikiの数式のとおりに実装できてはいるはずです。合っているのかはよくわかりません。

2016/9/14 追記: 正規化係数を1/root(N)にする(ノーマライズする)とwolfram alphaの数値と一致するようです。下のプログラムは正規化係数がDFTは1、IDFTは1/Nになっています。要するにwolfram alphaもこちらも間違えてなかったようです。

ソースコード

wikiの数式をそのままプログラムにしただけなので、説明はwiki読んでください。

#include <iostream>
#include <vector>
#include <map>
#include <cmath>
#include <complex>

using namespace std;
using comp = complex<double>;

vector<comp> fourier(const vector<comp>& f);
vector<comp> inverse_fourier(const vector<comp>& f);

int main() {
    vector<comp> f = {
        comp(0, 0), comp(1, 0), comp(4, 0), comp(9, 0), comp(16, 0),
        comp(25, 0), comp(36, 0), comp(49, 0), comp(64, 0), comp(81, 0),
        comp(100, 0), comp(121, 0), comp(144, 0), comp(169, 0), comp(196, 0),
        comp(225, 0), comp(256, 0), comp(289, 0), comp(324, 0), comp(361, 0),
        comp(400, 0), comp(441, 0), comp(484, 0), comp(529, 0), comp(576, 0),
        comp(625, 0), comp(676, 0), comp(729, 0), comp(784, 0), comp(841, 0),
        comp(900, 0), comp(961, 0), comp(1024, 0), comp(1089, 0), comp(1156, 0),
        comp(1225, 0), comp(1296, 0), comp(1369, 0)
    };
    vector<comp> res = fourier(f);

    vector<comp>::iterator itr;
    for (itr = res.begin(); itr != res.end(); ++itr) {
        cout << *itr << endl;
    }
    return 0;
}
vector<comp> fourier(const vector<comp>& f) {
    const int N = f.size();
    const double PI = 3.14159265358979265358979;
    const std::complex<double> i(0, 1);
    int t, x;
    vector<comp> F(N);
    for (t = 0; t < N; ++t) {
        F[t] = 0;
        for (x = 0; x < N; ++x) {
            F[t] += f[x] * exp(-i * (2.0 * PI * t * x / N));
        }
    }
    return F;
}
vector<comp> inverse_fourier(const vector<comp>& F) {
    const int N = F.size();
    const double PI = 3.14159265358979265358979;
    const std::complex<double> i(0, 1);
    vector<comp> f(N);
    int x, t;

    for (x = 0; x < N; ++x) {
        f[x] = 0;
        for (t = 0; t < N; ++t) {
            f[x] += F[t] * exp(i * (2 * PI * x * t / N));
        }
        f[x] /= N;
    }
    return f;
}

FFTも実装しました

ついでにFFTも実装しました。これと上のプログラムは出力が同じです。wolfram alphaとはずれています。

2016/9/14 追記: こちらも合っています

こっちについても、wikiの数式をそのままプログラムにしただけなので解説はしません。

ソースコードはこんな感じです

vcomp fourier(const vcomp& f) {
    int N = f.size();
    int P, Q, p, q, r, s;
    int i;
    const double PI = 3.14159265358979265358979;

    P = 1; Q = N;
    for (i = 2; i * i <= N; ++i) {
        if (N % i == 0) {
            P = i;
            Q = N / i;
        }
    }

    // step1
    vector<vcomp> f1(Q, vcomp(P));
    for (p = 0; p < P; ++p) {
        for (r = 0; r < Q; ++r) {
            f1[r][p] = 0;
            for (q = 0; q < Q; ++q) {
                const double theta = -2.0 * PI * r * q / Q;
                f1[r][p] += f[q * P + p] * comp(cos(theta), sin(theta));
            }
            const double theta = -2.0 * PI * p * r / P / Q;
            f1[r][p] *= comp(cos(theta), sin(theta));
        }
    }

    // step2
    vcomp F(N);
    if (P == 1) {
        for (r = 0; r < Q; ++r) {
            // s === 0
            F[r] = f1[r][0];
        }
    } else {
        for (r = 0; r < Q; ++r) {
            vcomp res = fourier(f1[r]);
            for (s = 0; s < P; ++s) {
                F[s * Q + r] = res[s];
            }
        }
    }
    return F;
}
vcomp inverse_fourier(const vcomp& f) {
    const int N = f.size();
    vcomp conj_f;

    vcomp::const_iterator const_itr;
    vcomp::iterator itr;
    for (const_itr = f.begin(); const_itr != f.end(); ++const_itr) {
        conj_f.push_back(conj(*const_itr));
    }

    vcomp res = fourier(conj_f);
    for (itr = res.begin(); itr != res.end(); ++itr) {
        *itr = conj(*itr);
        *itr /= N;
    }
    return res;
}