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

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

新規登録して質問してみよう
ただいま回答率
85.50%
C++

C++はC言語をもとにしてつくられた最もよく使われるマルチパラダイムプログラミング言語の1つです。オブジェクト指向、ジェネリック、命令型など広く対応しており、多目的に使用されています。

Q&A

解決済

3回答

5320閲覧

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

pandanoir

総合スコア72

C++

C++はC言語をもとにしてつくられた最もよく使われるマルチパラダイムプログラミング言語の1つです。オブジェクト指向、ジェネリック、命令型など広く対応しており、多目的に使用されています。

3グッド

0クリップ

投稿2016/09/09 02:07

編集2016/09/13 14:10

高速フーリエ変換のプログラムを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が正しくて私が間違えているのかもしれません。

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

cpp

1#include <iostream> 2#include <vector> 3#include <map> 4#include <cmath> 5#include <complex> 6 7using namespace std; 8using comp = complex<double>; 9using vcomp = vector<comp>; 10 11vcomp fourier(const vcomp& f); 12vcomp inverse_fourier(const vcomp& f); 13 14int main() { 15 vcomp f = { 16 comp(0, 0), comp(1, 0), comp(4, 0), comp(9, 0), comp(16, 0), 17 comp(25, 0), comp(36, 0), comp(49, 0), comp(64, 0) 18 }; 19 vcomp res = fourier(f); 20 21 vcomp::iterator itr; 22 for (itr = res.begin(); itr != res.end(); ++itr) { 23 cout << itr->real() << ", " << itr->imag() << endl; 24 } 25 return 0; 26} 27vcomp fourier(const vcomp& f) { 28 int N = f.size(); 29 int P, Q, p, q, r, s; 30 int i; 31 const double PI = 3.14159265358979265358979; 32 33 P = 1; Q = N; 34 for (i = 2; i * i <= N; ++i) { 35 if (N % i == 0) { 36 P = i; 37 Q = N / i; 38 } 39 } 40 41 // step1 42 vector<vcomp> f1(P, vcomp(Q)); 43 for (p = 0; p < P; ++p) { 44 for (r = 0; r < Q; ++r) { 45 f1[p][r] = 0; 46 if (Q == 2) { 47 f1[p][r] += f[p] * comp(1, 0); 48 f1[p][r] += f[P + p] * comp(r == 0 ? 1 : -1, 0); 49 } else if (Q == 4) { 50 const vcomp _e = {comp(1, 0), comp(0, -1), comp(-1, 0), comp(0, 1)}; 51 for (q = 0; q < Q; ++q) { 52 const int x = r * q % 4; 53 f1[p][r] += f[q * P + p] * _e[x]; 54 } 55 } else { 56 for (q = 0; q < Q; ++q) { 57 const double theta = -2 * PI * r * q / Q; 58 f1[p][r] += f[q * P + p] * comp(cos(theta), sin(theta)); 59 } 60 } 61 const double theta = -2 * PI * p * r / P / Q; 62 cout << p << " * " << r << " / " << P << " / " << Q << " = " << theta << endl; 63 f1[p][r] *= comp(cos(theta), sin(theta)); 64 } 65 } 66 67 // step2 68 vcomp F(N); 69 if (P == 1) { 70 for (r = 0; r < Q; ++r) { 71 // s === 0 72 F[r] = f1[0][r]; 73 } 74 } else { 75 vector<vcomp> f2(Q, vcomp(P)); 76 for (p = 0; p < P; ++p) { 77 for (r = 0; r < Q; ++r) { 78 f2[r][p] = f1[p][r]; 79 } 80 } 81 for (r = 0; r < Q; ++r) { 82 vcomp res = fourier(f2[r]); 83 for (s = 0; s < P; ++s) { 84 F[s * Q + r] = res[s]; 85 } 86 } 87 } 88 return F; 89} 90vcomp inverse_fourier(const vcomp& f) { 91 const int N = f.size(); 92 vcomp conj_f; 93 94 vcomp::const_iterator const_itr; 95 vcomp::iterator itr; 96 for (const_itr = f.begin(); const_itr != f.end(); ++const_itr) { 97 conj_f.push_back(conj(*const_itr)); 98 } 99 100 vcomp res = fourier(conj_f); 101 for (itr = res.begin(); itr != res.end(); ++itr) { 102 *itr = conj(*itr); 103 *itr /= N; 104 } 105 return res; 106}
Yoshikatsu_, mhashi, ikuwow👍を押しています

気になる質問をクリップする

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答3

0

ベストアンサー

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

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

python

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

投稿2016/09/14 02:25

pebble8888

総合スコア390

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

0

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

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

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

投稿2016/09/13 14:22

HogeAnimalLover

総合スコア4830

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

0

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

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

投稿2016/09/09 07:56

PineMatsu

総合スコア3579

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問