質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.58%

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

解決済

回答 3

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 3,360

pandanoir

score 65

高速フーリエ変換のプログラムを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;
}
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 3

checkベストアンサー

0

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])
>>> 

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

0

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

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

0

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

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

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.58%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る