高速フーリエ変換

高速フーリエ変換（こうそくフーリエへんかん、英: fast Fourier transform, FFT）は、離散フーリエ変換（英: discrete Fourier transform, DFT）を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換（英: inverse fast Fourier transform, IFFT）と呼ぶ。

歴史

高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した とされているクーリー–テューキー型FFTアルゴリズムを呼ぶ。同時期に高橋秀俊がクーリーとテューキーとは全く独立にフーリエ変換を高速で行うためのアルゴリズムを考案していた。しかし、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に対してどのように計算すればよいかについての説明をみることができる）。

クーリー–テューキー型FFTアルゴリズム

クーリー–テューキー型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。
分割統治法を使ったアルゴリズムで、N = N₁ N₂ のサイズの変換を、より小さいサイズである N₁,  N₂ のサイズの変換に分割していくことで高速化を図っている。
最もよく知られたクーリー–テューキー型アルゴリズムは、ステップごとに変換のサイズをサイズ N/2 の2つの変換に分割するので、2 の累乗次数に限定される。しかし、一般的には次数は 2 の累乗にはならないので、素因数が偶数と奇数とで別々のアルゴリズムに分岐する。
伝統的なFFTの処理実装の多くは、再帰的な処理を、系統だった再帰をしないアルゴリズムにより実現している。
クーリー–テューキー型アルゴリズムは変換をより小さい変換に分解していくので、後述のような他の離散フーリエ係数のアルゴリズムと任意に組み合わせることができる。とりわけ、N ≤ 8 あたりまで分解すると、固定次数の高速なアルゴリズムに切り替えることが多い。

プログラムの例

以下は、高速フーリエ変換のプログラムを Q = 4 の場合にMicrosoft Visual Basicの文法を用いて書いた例である。
Const pi As Double = 3.14159265358979   '円周率Dim Ndeg As Long '4^degDim 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^degReDim 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 j0Next 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 の間の部分)の繰り返し回数が Ndeg log₄ Ndeg となっている。