高速フーリエ変換のプログラムの計算結果がおかしい

Question

高速フーリエ変換のプログラムをwikiの数式に従って書いてみました。が、Wolfram Alphaの計算結果とずれてしまいます。具体的には、実部が3倍、虚部が-3倍ずれています。何がおかしいのでしょうか?ちなみに、Wolfram AlphaにはFourier[{0, 1, 4, 9, 16, 25, 36, 49, 64}]と打ち込みました。

行った検証

なんか文句つけられたので足掻いたあと記しときますね

まずwikiの数式を読み、解説を読みました。で、その通りに実装してみました。以下のコード中でstep1、step2とあるところがwikiと対応しています。また、離散フーリエ変換の定義式を元にプログラムを組み、それと比較する実験も行いました。そのプログラムと以下のプログラムは結果が一致しました。

次に、数式を小さいケースで手計算してみました。
Fourier[{1, 4, 9, 16}] というケースです。このフーリエ変換で得られるF(0)は、離散フーリエ変換の定義よりf(0)+f(1)+f(2)+f(3)つまり30となります。しかし、wolfram alphaでは15となっていました。どうもwolfram alphaが間違えている気がします。

しかし、私は別に数学者でもなんでもない一個人です。wolfram alphaが正しくて私が間違えているのかもしれません。

だから、プログラムの検証を求めます。これでもまだダメでしょうか？バカで幼稚だからよく分かりません

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

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

vcomp fourier(const vcomp& f);
vcomp inverse_fourier(const vcomp& f);

int main() {
    vcomp 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)
    };
    vcomp res = fourier(f);

    vcomp::iterator itr;
    for (itr = res.begin(); itr != res.end(); ++itr) {
        cout << itr->real() << ", " << itr->imag() << endl;
    }
    return 0;
}
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(P, vcomp(Q));
    for (p = 0; p < P; ++p) {
        for (r = 0; r < Q; ++r) {
            f1[p][r] = 0;
            if (Q == 2) {
                f1[p][r] += f[p] * comp(1, 0);
                f1[p][r] += f[P + p] * comp(r == 0 ? 1 : -1, 0);
            } else if (Q == 4) {
                const vcomp _e = {comp(1, 0), comp(0, -1), comp(-1, 0), comp(0, 1)};
                for (q = 0; q < Q; ++q) {
                    const int x = r * q % 4;
                    f1[p][r] += f[q * P + p] * _e[x];
                }
            } else {
                for (q = 0; q < Q; ++q) {
                    const double theta = -2 * PI * r * q / Q;
                    f1[p][r] += f[q * P + p] * comp(cos(theta), sin(theta));
                }
            }
            const double theta = -2 * PI * p * r / P / Q;
            cout << p << " * " << r << " / " << P << " / " << Q << " = " << theta << endl;
            f1[p][r] *= comp(cos(theta), sin(theta));
        }
    }

    // step2
    vcomp F(N);
    if (P == 1) {
        for (r = 0; r < Q; ++r) {
            // s === 0
            F[r] = f1[0][r];
        }
    } else {
        vector<vcomp> f2(Q, vcomp(P));
        for (p = 0; p < P; ++p) {
            for (r = 0; r < Q; ++r) {
                f2[r][p] = f1[p][r];
            }
        }
        for (r = 0; r < Q; ++r) {
            vcomp res = fourier(f2[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;
}

Answer 1

3倍なので、normalizationがかかっているだけじゃないでしょうか。要素数が9なので√9=3ですし。
normalizationについてはここらに記載がありました。
http://docs.scipy.org/doc/numpy/reference/routines.fft.html#module-numpy.fft

pythonのnumpy.fftで計算したらこうなりました。

$ python
>>> import numpy
>>> numpy.fft.fft([0,1,4,9,16,25,36,49,64])
array([ 204.00000000  +0.j        ,   -2.03115523+111.27283549j,
        -29.60875519 +48.2660205j ,  -34.50000000 +23.3826859j ,
        -35.86008958  +7.14124272j,  -35.86008958  -7.14124272j,
        -34.50000000 -23.3826859j ,  -29.60875519 -48.2660205j ,
         -2.03115523-111.27283549j])
>>> 

Answer 2

数値演算は誤差を避けられないため、単純な検証は難しいです。とりあえず高速でない離散フーリエ（DFT）を先に作り、これと比較するという検証方法が考えられます。

ちなみに、複素数計算とDFTを先に定義しておけば、高速フーリエ（FFT）は再帰関数を用いて非常に単純に記述することができます。

ヒント
要素数が奇数である場合：DFTに丸投げ
要素数が偶数である場合：
１ 奇数番目のみを再帰的にFFT→この結果をFFT１と呼ぶ
２ 偶数番目のみを再帰的にFFT→この結果をFFT２と呼ぶ
３ FFT１とFFT２、回転子を用いてFFTの結果となる配列を作る

Answer 3

一般的に、数値演算の検証は難しいです。
ある条件ではうまくいけても条件が変わればおかしいなどかなりの検証が必要です。できれば、信頼できる検証済みのFFTライブラリの結果と突き合わせて検証するなどが必要です。

まずは「C++ FFT ソース」でググればいくつか出て来るので、それらのソースと見比べてみればどうでしょうか？