出典(authority):フリー百科事典『ウィキペディア(Wikipedia)』「2014/05/28 16:54:53」(JST)
高速フーリエ変換(こうそくフーリエへんかん、英: Fast Fourier Transform、FFT)とは、離散フーリエ変換 (Discrete Fourier Transform、DFT) を計算機上で高速に計算するアルゴリズム。逆変換をIFFT (Inverse FFT) という。
高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー(英語版) (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した[1]とされているCooley-Tukey型FFTアルゴリズム(英語版)を呼ぶ[2]。しかし、1805年ごろにガウスが同様のアルゴリズムを独自に発見していた[3]。
複素関数f(x)の離散フーリエ変換 (DFT)である複素関数F(t) は以下で定義される。
このとき、を標本点と言う。
これを直接計算したときの時間計算量は O (N2 ) である(Oはランダウの記号)。
高速フーリエ変換は、この結果を、次数 N が2の累乗のときに O (N log N ) の計算量で得るアルゴリズムである。より一般的には、次数が と素因数分解できるとき、 の計算量となる。次数が2の累乗のときが最も高速に計算でき、アルゴリズムも単純になるので、0詰めで次数を調整することもある。
高速フーリエ変換を使って、畳み込み積分などの計算を高速に求めることができる。これも計算量をO (N2 ) → O (N log N ) に落とせる。
現在は、初期の手法[1]をより高速化したアルゴリズムが使用されている。
逆変換 (IFFT) は正変換 (FFT) と同じと考えて良いが、指数の符号が逆であり、係数が掛かる。
離散フーリエ変換を以下で定義したとき、
逆変換は、
となる。
このため、の離散フーリエ逆変換を求めるには、
とすれば良く、正変換の高速フーリエ変換のプログラムがあれば、逆変換は容易に作ることができる。
詳細は「Cooley-Tukey型FFTアルゴリズム(英語版)」を参照
Cooley-Tukey型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。
分割統治法 (divide and conquer) を使ったアルゴリズムで、 N = N1 N2 のサイズの変換を、より小さいサイズである N1 、N2 のサイズの変換に分割していくことで高速化を図っている。
最もよく知られたCooley-Tukey型アルゴリズムは、ステップごとに変換のサイズをサイズ N / 2 の2つの変換に分割するので、2の累乗次数に限定される。しかし、一般的には次数は2の累乗にはならないので、素因数が偶数と奇数とで別々のアルゴリズムに分岐する。
伝統的なFFTの処理実装の多くは、再帰的な処理を、系統だった再帰をしないアルゴリズムにより実現している。
Cooley-Tukey型アルゴリズムは変換をより小さい変換に分解していくので、後述のような他の離散フーリエ係数のアルゴリズムと任意に組み合わせることができる。とりわけ、N ≤ 8 あたりまで分解すると、固定次数の高速なアルゴリズムに切り替えることが多い。
離散フーリエ係数は、1の原始 N 乗根の1つ を使うと、次のように表せる。
例えば、N = 4 のとき、離散フーリエ係数は行列を用いて表現すると(W ≡ W4 と略記)
となる。入力列 xk を添字の偶奇で分けて、以下のように変形する。
すると、サイズ2のFFTの演算結果を用いて表現でき、サイズの分割ができる。
また、この分割手順を図にすると蝶のような図になることから、バタフライ演算とも呼ばれる。
バタフライ演算は、計算機上ではビット反転で実現される。DSPの中には、このバタフライ演算のプログラムを容易にするため、ビット反転アドレッシングを備えているものがある。
N =PQ とする。N 次離散フーリエ変換:
を、n = 0, 1, ... , N -1 について計算することを考える。n , k を次のように書き換える。ただし0 ≤ n ≤ N -1、0 ≤ k ≤ N -1である。
すると
ここで、
と置くと、
となる。即ち、の計算は、次の2ステップになる。
ステップ1、2は、N = PQ 次の離散フーリエ変換を、Q 次の離散フーリエ変換と回転因子の掛け算の実行により、Q 組(r = 0, 1, ... , Q -1)のP 次離散フーリエ変換に分解したと見ることができる。
N = Qk (P = Qk -1 )の場合には、上を繰り返せば、Q 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換(何もしないことと同じ)にまで下げると、F (t )を求めることができる。特に、Q が2または4の場合は、Q 次の離散フーリエ変換は非常に簡単な計算になる。
このため、2の累乗あるいは4の累乗次の離散フーリエ変換は簡単に計算できる。実務的に用いられるのは、Q = 2かQ = 4の場合のみである。なお、Q = 2かQ = 4の場合のこの部分のQ 次の離散フーリエ変換のことを、バタフライ演算と言う。
また、Q = 2かQ = 4の場合において、計算を終了するまでに何回の「掛け算」が必要かを考える。符号の逆転、実部虚部の交換は「掛け算」として数えなければ、回転因子の掛け算のみが「掛け算」である。N = Qkの次数を1落とすためにN回の「掛け算」が必要であり、次数をkから0に落とすにはそれをk回繰り返す必要があるため、「掛け算」の数はNk=NlogQNとなる。高速フーリエ変換の計算において時間がかかるのは「掛け算」の部分であるため、これが「高速フーリエ変換では計算速度は O (NlogN ) になる」という根拠になっている。
上記の説明で、N = Qk (P = Qk -1 )の場合、N = Qk 個のデータから、N = Qk 個の計算結果
を計算する場合に、メモリの節約のため、0≤q ≤Q -1、0≤r ≤Q -1 を利用し、計算結果f1(p ,r ) を元データ のあった場所に格納することが多い。これが次の次数Qk -1 でも繰り返されるため、とすると、次の次数の計算結果はのあった場所に格納される。繰り返せば、とすると、計算結果はのあった場所に格納される。
一方、
を、r を固定しs を変数としたQk -1 次離散フーリエ変換と見なして、とすると、
となる。繰り替えせば、
となるが、左辺について
より、また右辺について
より。このため、
これはのあった場所に格納されている。
このように、求める解 が のあった場所に格納されていることを、ビット反転と言う。これは、Q 進数で表示した場合、 は となるのに対し、は逆から読んだとなるためである。
Microsoft Visual Basicの文法で、高速フーリエ変換のプログラムを書くと、次のような例が考えられる。この例ではQ = 4である。
(プログラムはじまり) Const pi As Double = 3.14159265358979 '円周率 Dim Ndeg As Long '4^deg Dim Pdeg As Long '4^(deg-i) Dim CR() As Double '入力実数部 Dim CI() As Double '入力虚数部 Dim FR() As Double '出力実数部 Dim FI() As Double '出力虚数部 deg=5 '任意に設定。5ならN=4^5=1024で計算 Ndeg=4^deg ReDim CR(Ndeg - 1) As Double '入力実数部 ReDim CI(Ndeg - 1) As Double '入力虚数部 ReDim FR(Ndeg - 1) As Double '出力実数部 ReDim FI(Ndeg - 1) As Double '出力虚数部 'ここで、変換される関数の実部をCR(0)からCR(Ndeg-1)に、虚部をCI(0)からCI(Ndeg-1)に入力しておくこと 'フーリエ変換 For i = 1 To deg Pdeg = 4 ^ (deg - i) For j0 = 0 To 4 ^ (i - 1) - 1 For j1 = 0 To Pdeg - 1 j = j1 + j0 * Pdeg * 4 'バタフライ演算(Q=4) w1 = CR(j) + CR(j + Pdeg) + CR(j + 2 * Pdeg) + CR(j + 3 * Pdeg) w2 = CI(j) + CI(j + Pdeg) + CI(j + 2 * Pdeg) + CI(j + 3 * Pdeg) w3 = CR(j) + CI(j + Pdeg) - CR(j + 2 * Pdeg) - CI(j + 3 * Pdeg) w4 = CI(j) - CR(j + Pdeg) - CI(j + 2 * Pdeg) + CR(j + 3 * Pdeg) w5 = CR(j) - CR(j + Pdeg) + CR(j + 2 * Pdeg) - CR(j + 3 * Pdeg) w6 = CI(j) - CI(j + Pdeg) + CI(j + 2 * Pdeg) - CI(j + 3 * Pdeg) w7 = CR(j) - CI(j + Pdeg) - CR(j + 2 * Pdeg) + CI(j + 3 * Pdeg) w8 = CI(j) + CR(j + Pdeg) - CI(j + 2 * Pdeg) - CR(j + 3 * Pdeg) CR(j) = w1 CI(j) = w2 CR(j + Pdeg) = w3 CI(j + Pdeg) = w4 CR(j + 2 * Pdeg) = w5 CI(j + 2 * Pdeg) = w6 CR(j + 3 * Pdeg) = w7 CI(j + 3 * Pdeg) = w8 '回転因子 For k = 0 To 3 w1 = Cos(2 * pi * j * k / Pdeg / 4) w2 = -Sin(2 * pi * j * k / Pdeg / 4) w3 = CR(j + k * Pdeg) * w1 - CI(j + k * Pdeg) * w2 w4 = CR(j + k * Pdeg) * w2 + CI(j + k * Pdeg) * w1 CR(j + k * Pdeg) = w3 CI(j + k * Pdeg) = w4 Next k Next j1 Next j0 Next i 'ビット反転 For i = 0 To Ndeg - 1 k = i k1 = 0 For j = 1 To deg k1 = k1 + (k - Int(k / 4) * 4) * 4 ^ (deg - j) k = Int(k / 4) Next j FR(i) = CR(k1) FI(i)=CI(k1) Next i (プログラムおわり)
この例では、最深部(For k、Next kの間の部分)の繰り返し回数がNdegNdegとなっている。
この節の加筆が望まれています。 |
多くの実用面において、FFTへの入力は実数のみの列(実入力)であり、このとき出力の列は次の対称性を満たす(*は複素共役):
そして、多くの効率的なFFTアルゴリズム(例えば[4])はこの条件を前堤に設計されている。
実入力への効率化には、次のような手段がある。
かつては実入力に対するフーリエ係数は離散ハートリー変換(英語版) (DHT : Discrete Hartley transform) でより効率的に計算できると思われていた。しかし、その後最適化された離散フーリエ変換(DFT:Discrete Fourier Transform)アルゴリズムが、離散ハートリー変換アルゴリズムよりも必要とする演算回数を下まわることがわかった。また、Bruun's FFTアルゴリズムは、はじめこそ実入力に対して有利といわれていたが、今ではそういわれてはいない。
さらに偶奇の対称性を持つ実入力の場合にはさらに最適化ができ、DFTはDCTやDST(英語版)になり、時間とメモリーの(ほとんど)2つに関して最適化が得られる。この場合には直接DFTのアルゴリズムを応用しないでも、DCTやDSTの計算によってフーリエ係数を求めることができる。
[ヘルプ] |
全文を閲覧するには購読必要です。 To read the full text you will need to subscribe.
リンク元 | 「FFT」「fast Fourier transform」「fast Fourier transformation」 |
関連記事 | 「変換」「高速」「フーリエ変換」「フーリエ」 |
.