高速フーリエ変換

高速フーリエ変換（こうそくフーリエへんかん、）は、離散フーリエ変換（）を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換（）と呼ぶ。
概要.
複素関数 の離散フーリエ変換である複素関数 は以下で定義される。
このとき、 を標本点と言う。
これを直接計算したときの時間計算量は、ランダウの記号を用いて表現すると である。
高速フーリエ変換は、この結果を、次数が2の累乗のときに の計算量で得るアルゴリズムである。より一般的には、次数が と素因数分解できるとき、 の計算量となる。次数が の累乗のときが最も高速に計算でき、アルゴリズムも単純になるので、 詰めで次数を調整することもある。
高速フーリエ変換を使って、畳み込み積分などの計算を高速に求めることができる。これも計算量を から まで落とせる。
現在は、初期の手法 をより高速化したアルゴリズムが使用されている。
逆変換.
逆変換は正変換と同じと考えて良いが、指数の符号が逆であり、係数 が掛かる。
高速フーリエ変換のプログラム中、どの符号が逆転するかを一々分岐させると、分岐の判定に時間がかかり、パフォーマンスが落ちる。一方、正変換のプログラムと、逆変換のプログラムを両方用意しておくことも考えられるが、共通部分が多いため、無駄が多くなる。このため、複素共役を使った次のような方法が考えられる。
離散フーリエ変換を
で定義したとき、逆変換は
となる。
このため、 の離散フーリエ逆変換を求めるには、
とすれば良く、正変換の高速フーリエ変換のプログラムがあれば、逆変換は容易に作ることができる。
アルゴリズム.
クーリー–テューキー型FFTアルゴリズム.
クーリー–テューキー型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。
分割統治法を使ったアルゴリズムで、 のサイズの変換を、より小さいサイズである , のサイズの変換に分割していくことで高速化を図っている。
最もよく知られたクーリー–テューキー型アルゴリズムは、ステップごとに変換のサイズをサイズ の2つの変換に分割するので、 の累乗次数に限定される。しかし、一般的には次数は の累乗にはならないので、素因数が偶数と奇数とで別々のアルゴリズムに分岐する。
伝統的なFFTの処理実装の多くは、再帰的な処理を、系統だった再帰をしないアルゴリズムにより実現している。
クーリー–テューキー型アルゴリズムは変換をより小さい変換に分解していくので、後述のような他の離散フーリエ係数のアルゴリズムと任意に組み合わせることができる。とりわけ、 あたりまで分解すると、固定次数の高速なアルゴリズムに切り替えることが多い。
原理の簡単な説明.
離散フーリエ係数は、の原始 乗根の1つ を使うと、次のように表せる。
例えば、 のとき、formula_5、formula_6とすれば、離散フーリエ係数は行列を用いて表現すると（ と略記）
となる。入力列 を添字の偶奇で分けて、以下のように変形する。
すると、サイズ のFFTの演算結果を用いて表現でき、サイズの分割ができる。
また、この分割手順を図にすると蝶のような図になることから、バタフライ演算とも呼ばれる。
バタフライ演算は、計算機上ではビット反転で実現される。DSPの中には、このバタフライ演算のプログラムを容易にするため、ビット反転アドレッシングを備えているものがある。
原理の説明.
 とする。 次離散フーリエ変換を、以下のように 次離散フーリエ変換と 次離散フーリエ変換に分解する。
 次離散フーリエ変換：
を、 について計算することを考える。, を次のように書き換える。ただし また である。
すると
ここで、
と置くと、
となる。即ち、 の計算は、次の2ステップになる。
ステップ1、2は、 次の離散フーリエ変換を、 次の離散フーリエ変換と回転因子の掛け算の実行により、 組 () の 次離散フーリエ変換に分解したと見ることができる。
 () の場合には、上を繰り返せば、 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換（何もしないことと同じ）にまで下げると、を求めることができる。特に、 がまたはの場合は、次の離散フーリエ変換は非常に簡単な計算になる。
このため、の累乗あるいはの累乗次の離散フーリエ変換は簡単に計算できる。実務的に用いられるのは、 か の場合のみである。なお、 か の場合のこの部分の次の離散フーリエ変換のことを、バタフライ演算と言う。
また、かの場合において、計算を終了するまでに何回の「掛け算」が必要かを考える。符号の逆転、実部虚部の交換は「掛け算」として数えなければ、回転因子の掛け算のみが「掛け算」である。の次数を1落とすために回の「掛け算」が必要であり、次数をからに落とすにはそれを回繰り返す必要があるため、「掛け算」の数は となる。高速フーリエ変換の計算において時間がかかるのは「掛け算」の部分であるため、これが「高速フーリエ変換では計算速度は になる」ことの根拠になっている。
ビットの反転.
上記の説明で、formula_19 の場合、 個のデータformula_20から、 個の計算結果
を計算する場合に、メモリの節約のため、 と を利用し、計算結果 formula_22 を元データformula_23 のあった場所に格納することが多い。これが次の次数 でも繰り返されるため、formula_24とすると、次の次数の計算結果formula_25はformula_26のあった場所に格納される。繰り返せば、formula_27とすると、計算結果formula_28はformula_29のあった場所に格納される。
一方、
を、 を固定し を変数とした 次離散フーリエ変換と見なして、formula_31とすると、
となる。繰り替えせば、
となるが、左辺について
より , また右辺について
より 。このため、
これは formula_37 のあった場所に格納されている。
このように、求める解 formula_38 が formula_37 のあった場所に格納されていることを、ビット反転と言う。これは、 進法で表示した場合、formula_40 は formula_41となるのに対し、formula_42は逆から読んだformula_43となるためである。
プログラムの例.
以下は、高速フーリエ変換のプログラムを の場合にMicrosoft Visual Basicの文法を用いて書いた例である。
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
この例では、最深部 (codice_1、codice_2 の間の部分)の繰り返し回数が codice_3 codice_3 となっている。
実数および対称的な入力への最適化.
多くの応用において、FFTに対する入力データは実数の列(実入力)であり、このとき変換された出力の列は次の対称性を満たす（ は複素共役）：
そこで、多くの効率的なFFTアルゴリズム は入力データが実数であることを前提に設計されている。
入力データが実数の場合の効率化の手段としては、次のようなものがある。
かつては実数の入力データに対するフーリエ係数を求めるのには、実数計算だけで行える (discrete Hartley transform, DHT)を用いると効率的であろうと思われていた。しかしその後に、最適化された離散フーリエ変換 (discrete Fourier transform, DFT) アルゴリズムの方が、離散ハートリー変換アルゴリズムに比べて必要な演算回数が少ないということが判明した。また当初は、実数入力に対してブルーン (Bruun) FFT アルゴリズムは有利であると云われていたが、その後そうではないことが判った。
また、偶奇の対称性を持つ実入力の場合には、DFTはDCTやとなるので、演算と記憶に関してほぼ2倍の効率化が得られる。よって、そのような場合にはDFTのアルゴリズムをそのまま適用するよりも、DCTやDSTを適用してフーリエ係数を求める方が効率的である。
歴史.
高速フーリエ変換といえば一般的には1965年、 (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した とされているを呼ぶ。同時期に高橋秀俊がクーリーとテューキーとは全く独立にフーリエ変換を高速で行うためのアルゴリズムを考案していた。しかし、1805年頃に既にガウスが同様のアルゴリズムを独自に発見していた（本ページの外部リンク先に同じ文章PDFへのリンクがある）。ガウスの論文以降、地球物理学や気候や潮位解析などの分野などで測定値に対する調和解析は行われていたので、計算上の工夫を必要とする応用分野で受け継がれていたようである（たとえば、Robart L. Nowack: "Development of the FFT and Applications in Geophysics", in Proceedings of the Cornelius Lanczos International Centenary Conference,SIAM, ISBN 978-0898713398 (1994), pp.395--397、の中では Danielson and Lanczos（1942年）などの先行例をあげている。和書でも沼倉三郎：「測定値計算法」、森北出版、（1956年）には，一般の合成数Nに対してではないが人が計算を行う場合にある程度の大きさの合成数Nに対してどのように計算すればよいかについての説明をみることができる）。
以下の書籍にも、天体観測の軌道の補間のためにガウスが高速フーリエ変換を利用したことが書かれている。
学習用図書.
今後記述を追加の予定