FFTプログラム（1次元）を使って高次元のFFTを計算するプログラム

以下、長文になりますがご容赦。
　1次元のFFTのフォートランソースプログラムはネットにも本の付録にもあります。高次元のソースプログラム（2次元、3次元）を作りたいと思っています。私は以前にもこの質問をしており、1次元のプログラムを何度も繰り返して作成できる、ということが分かりました。
　添付した画像はこのような方向かなと思って式を調べたものです。画像のようにｘ方向への変換（赤字）とその後ｙ方向への変換（そのあとの黒字式3行目）（都合2回）ということだろうと思います。これによりますと確かにできることがわかります。（このような理解でもいいでしょうか？）
　それでも少し疑問が残ります。赤字の変換は変換対象 f(x,y) が実数であり、その後の黒字の変換はg(kx,y)を変換していることになりますが、こちらは複素数ですね。また、たとえば、X方向に64個のデータでFFTをかけると波数kxの数は32個で位相情報（sin, cos)が32個づつ出るから系列の個数は64にはなるということではないかと思います。
　ということなので、1次元の計算を2回する、ということですが、1回目と2回目で若干様子が違うということになります。ただし、f(x,y)は実数ですが、複素数の変数の実部（虚部はゼロ）に代入してFFTルーチンに入れたようでもあります。だんだんあやふやになります。またkx, kyの定義についても２π(パイ)をかけるものがあったり、出てきたフーリエ成分はデータ数Nが掛けられているとかプログラムの作り方が作成者に依存するところがあったりするので正しく動作するかどうか不明です。
　フーリエ変換のプログラムが正しく動作しているかどうかは、すぐさま逆変換して比較するわけですが、フーリエ変換と逆変換は言わば同じものであり、片方が正しく動作すれば逆もほぼ正しいでしょうし、間違い（変換）と間違い（逆変換）が相殺して正しく見せてしまっているかも知れません。
　ということで、私の希望としては、これで正しく動作するというセットがあると安心して使えるのに、ということです。ネットを探してみると、２次元はあるかもな、というところです。画像処理として利用されている例があるからですね。３次元はどうでしょうか。ある、と聞いていますが、見つかりません。もしあるということであれば、ダウンロード情報を教えて頂きたいのですが。
　ソースプログラムはフォートランだとありがたいですが、Cでも解読できると思います。
実験して試行錯誤して解決するというよりも、先達に指定して頂くとありがたいです。また、人それぞれ自分のプログラムを持つ必要がなく、動作確認がされているソースを１つ持ちたいと思っているのです。あるいは1次元のFFTを複数回使って間違いなく計算する方法を教えて戴きたいのですが。
長文で申し訳ありませんが、よろしくお願いします。

No.4 ベストアンサー 回答者： スチールウール 回答日時： 2015/10/02 15:53

No.3の回答へのお礼を受けて

＞まず、j=1に固定してi=1~64を1DのFFTにかけます。結果,i=1~64の成分が複素数で出力されます。それをg(i,1) (i=1~64)に保存し、j=2,3,4...64と行くと、g(i,j),(64x64)ができます(図の赤字部)。次に、i=1に固定してg(1,j),j=1~64を1DのFFTにかけ、s(1,j),j=1~64に格納し、i=2,3,4..64と進めると、s(i,j) (64x64)が算出され、２次元フーリエ変換となる、という考え方です。
はい、この考え方で合っています

ソースコードは多分、前半部分がbit逆順配列への並べ替え、後半部分がFFTだと思います
コメントがないので解りづらいですが、頭のなかで計算する限りあっていると思います

この回答へのお礼

回答有難うございます。
また、長くお付き合い頂き、有難うございました。２次元に関しては、勝手に作った２次元のデータを上記のように処理してフーリエ成分を計算した後、その逆変換を２回やったら元の勝手に作ったデータが出てきたのでうまく行っているんだろうなあとは思います。３次元も実施します。
　実は最終目的は、物理空間の数値計算（時間発展）を波数空間のスペクトルの発展に置き換えて、波数空間の結果を逆変換して物理空間の様子を調べるというものです。スペクトル法という物理数値計算に属するのですが。

お礼日時：2015/10/10 00:23

No.3 回答者： スチールウール 回答日時： 2015/09/29 10:39

フーリエ変換と逆変換についてですが、基本はexp(-ikx)を掛けて積分するのでほぼ一緒です

何が変わってくるかというと、出てきた係数に最後に○倍するかどうかが変わります
昔習っただけなので余り詳しくは覚えていませんが、フーリエ変換でM倍と逆変換N倍して、MN=1/2πとならないと元信号に戻らないです
今回の添付の図ですと、すでにフーリエ変換で1/2π掛かっているので、逆変換ではこの項がいらないと思います
(そうでないなら、引数の1,-1も要らなくなるので)

＞その手持ちFFT（1次元）を繰り返し利用して2次元、3次元になるだろうか、というのが質問の主旨ですね。FFTのプログラムは1次元での動作が確認できたらそのことが、そのまま高次元に利用できる担保となるか、ということなのですが。
質問の意図がすみませんが汲み取れません
高次元FFT(例えば3D)のプログラムを組み、それに対して256x1x1のようなインプットを入れると言うことでしょうか？
それですと残念ですが担保？にはなっても保証にはなりません
一回目のフーリエ変換で返ってきた値(上の図でのg(kx,y))をきちんと次の入力にしているのかが分からないからです
また、fortranではなくCの場合、最高次元目はメモリ上に綺麗に並びますがその他の次元は並ばないため、ポインタ演算でstep値を設定する必要があります
一番いいのは、2次元三角関数 (sin(fx*x + fy*y) )や3次元三角関数を入れて結果を見ることです

1次元を転用したら簡単に高次元のプログラムを作れるかということでしたら、フーリエ変換の計算を省略できるので比較的簡単に出来ます
ただ、やはりどの次元のデータを1Dとして入力し、帰ってきた値をどこにどう保存するのか？などを新たに組む必要が有りそこで結構間違えます

＞FFTのプログラムはほぼ標準化されておりデータの格納の仕方まで含めてバラエティやクセがないということになるのでしょうか。
そんなことは無いです
wikipediaでFFTを検索してもらえばわかると思いますが、複数のアルゴリズムがあります
同じアルゴリズムでも、データの受け取り方(realを前半、imageを後半にするのか、別変数にするのか、最近のCだとcomplex型が使えるのでそれを使うのか)などで変わってきます
Cooley-Tukey型FFTアルゴリズムだと、逆向きで組んでしまうこともあります
Cooley-Tukeyの場合、2の冪乗のデータしか扱えないので後ろに0を挿入するのが一般的だと思いますが、前後にpaddingする人も居るでしょうし
FFTWの場合はたしか、3の倍数が混じっている(8x3=24など)場合も高速化するアルゴリズムが組み込まれているはずです

とりあえずは、多次元FFT(高速フーリエ変換)ではなく多次元FT(フーリエ変換)のプログラムを組まれてはどうでしょうか？

この回答へのお礼

誠に有難うございます。
２次元空間の実変数はf(i,j) でx=idx, y=jdyで、i,jの大=x,yの大です。まず、j=1に固定してi=1~64を1DのFFTにかけます。結果,i=1~64の成分が複素数で出力されます。それをg(i,1) (i=1~64)に保存し、j=2,3,4...64と行くと、g(i,j),(64x64)ができます(図の赤字部)。次に、i=1に固定してg(1,j),j=1~64を1DのFFTにかけ、s(1,j),j=1~64に格納し、i=2,3,4..64と進めると、s(i,j) (64x64)が算出され、２次元フーリエ変換となる、という考え方です。以下は1DのFFTです。考えに対応しているでしょうか。

subroutine fft(n,a,nm,iopt)
complex a(nm),b,c
j=1
do 4 i=1,n
if(i.ge.j) goto 1
b=a(j)
a(j)=a(i)
a(i)=b
1 m=n/2
2 if(j.le.m) goto 3
j=j-m
m=m/2
if(m.ge.2) goto 2
3 j=j+m
4 continue

nhold=1

5 if(nhold.ge.n) goto 8

int=2*nhold
do 7 k=1,nhold
p=3.141593*float(-iopt*(k-1))/float(nhold)
c=cmplx(0.0,p)
do 6 i=k,n,int
j=i+nhold
b=a(j)*cexp(c)
a(j)=a(i)-b
a(i)=a(i)+b
6 continue
7 continue
nhold=int
goto 5
8 if(iopt.eq.-1) goto 10
do 9 i=1,n
a(i)=a(i)/float(n)
9 continue
10 return
end

お礼日時：2015/10/02 12:58

No.2 回答者： スチールウール 回答日時： 2015/09/27 08:36

幾つかだけ回答を

まず、添付式の部分に関しての理解は正しいと思います
ただし、FFTは離散計算ですので本来ならフーリエ変換ではなくフーリエ級数展開で考えるべきですが

g(kx,y)が複素数になるとのことですが、こちらはそういうものだとしか言いようがないです
kx,kyに2πやNで割っていたりするのは実用的なことを考えているからです
2π/Nを使っていないと、例えば画像のサイズが512で合った場合、kx=0.1というと20π画素で一周期となり整数になりません
そこで、2π/Nする (どちらかというとkxにするのではなくxにしてる気がしますが)ことにより、画像内で何周期かという周波数が求まります
要するに、元のままだとkxの単位はcycles/pixとなり、2π/N掛けるとcyclex/Imageになり分かり易いのです

フーリエ変換と逆変換は同じとおっしゃられますが、正確には違います
今回の場合ですとフーリエ変換の方に1/(2π)が掛かっているため逆変換では必要ありません
逆変換と同じにしたい場合は係数を元々√にしておきます
また、どうしても計算機の場合誤差が出るため逆変換しても虚数部に少し値が残ることもありえます

FFTのソースコードに関しては、他人のものと比較してもあまり意味が無いと思います
なぜなら、FFTをする際に、データを2の冪乗にするためにzero paddingするのか、他の方法を使うのかなど幾つかアルゴリズムがあるからです
まだフーリエ変換のソースコードのほうが見つかりやすいと思います
基本的にFFT(Fast Fourier Transform：高速フーリエ変換)と言って高速化アルゴリズムを使っているのに、2Dや3Dに拡張した時にわざわざ書き換えているものは少ないと思います
あったとしてもインターフェース的なものであり内部では結局1D-FFTを呼んでいるだけということが多いです

この回答へのお礼

回答有り難うございます。今手持ちのFFTソースコードは本にサンプルコードとして載っているもので、実際に入力してみたら思ったように動作します（正も逆も）。
また、”フーリエ変換・逆変換が同じ”と私が言ったのは、そのサンプルコードにはスイッチ（－１とか１とか）があり、（少しは違いがありますが）同じようなループを回っているだけのように見える（FFTのソースは100行程度）ということです。
（話を戻しまして）
その手持ちFFT（1次元）を繰り返し利用して2次元、3次元になるだろうか、というのが質問の主旨ですね。FFTのプログラムは1次元での動作が確認できたらそのことが、そのまま高次元に利用できる担保となるか、ということなのですが。FFTのプログラムはほぼ標準化されておりデータの格納の仕方まで含めてバラエティやクセがないということになるのでしょうか。

お礼日時：2015/09/29 02:10

No.1 回答者： kmee 回答日時： 2015/09/23 08:49

https://oshiete.goo.ne.jp/qa/9027312.html
の続きですよね。

前にも言いましたけど、 64個の複素数をFFTにかけると、64個の複素数が出てきます。
実数を入力する場合は 虚数部が0の複素数として計算します。

スペクトル解析等では、半分の32個に注目する(残り半分は折り返しになっているので、データとして意味がない)という使い方はします。
ですが、FFT自体は64個の複素数を返します。

> X方向に64個のデータでFFTをかけると波数kxの数は32個で位相情報（sin, cos)が32個づつ出るから系列の個数は64にはなるということではないかと思います。

というのは間違っています。

ここでkxは波数ではなく、系列の個数です。

出力は複素数です。「位相情報（sin, cos)が32個づつ出るから系列の個数は64」等と分けて数えてはいけません。
実数の配列としてプログラムした場合は
・実数部用64、虚数部用64の計128個の実数を入力
・実数部用64、虚数部用64の計128個の実数が出力される
とするでしょう。ただし、どちらも系列の個数は「64」です。
この「プログラムでの実装方法における実数の数」と「系列の個数」を混同していませんか?

多次元で計算する場合は、1段目で求めた複素数をそのまま2段目で使います。


> ある、と聞いていますが、見つかりません。もしあるということであれば、ダウンロード情報を教えて頂きたいのですが。

どうやって探しました?
「3次元 FFT fortran」で検索したら
http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html
が一番上に出てきましたよ。
※ こちらで試したわけではないので、保証はしません。

この回答へのお礼

詳細な回答有難うございます。とにかく逆も正も複素数を入力して複素数を返すということですね。
実数の場合は虚部をゼロにするということで。
以下確認ですが。
n=64でx方向にΔxで分解された空間の広がり(nΔx)のデータがあるとします。入力としては、64個の複素数（実部64個（実空間のデータ）、虚部64個で値はゼロ）。
スペクトルの分解能は2π/(nΔx)=Δkであり、1Δk, 2Δk, ３Δk,....64Δkに対応した複素数が出力されるということですね（これも実部64個、虚部64個）。複素数は実部がｘで、虚部がiyであり、x^2+y^2がそれぞれの波数に対応したパワー、atan(y/x)が位相θになる。ちょっと気になるのが定数成分ですが。虚部がゼロの定数成分と63個の複素数ということでしょうか。これをそっくりそのままもう１回y方向にフーリエ変換すればそれでよし、ということで。作り方に個性が出そうな気もしますが。

プログラムは作りこみ方（実装、インプリメント）が人によって異なるようなことが起こるので、複数の人が作ると理屈が分っていたとしてもミスが生じる可能性があるわけで、すべてを通しで作って（かつサブルーチン化して）もらうのが一番ありがたいです。サイトの方にもアクセスしてみます。

お礼日時：2015/09/23 11:14