高速フーリエ変換とは (コウソクフーリエヘンカンとは) [単語記事]

高速フーリエ変換(Fa st Fourier Transform、FFT)とは、高速な離散フーリエ変換である。

概要

高速フーリエ変換とは、離散フーリエ変換(DFT)を高速に行うためのアルゴリズムの総称である。DFTはO(N²)程度の計算量が要求されるが、FFTならばO(Nl ogN)程度まで計算量を削減できる。どれくらいかと言うと、通常のDFTは4096個のデータで100秒くらいかかるが、FFTなら20ミリ秒くらいで終わってしまうほど。

プログラムで実装する場合はクイックソートのような再帰的プログラムが簡潔になる。具体的には、データの列から半分の長さの列を2つ作り、それぞれの列に再帰的にプログラムを適用させる。

以下に記載する基本原理ではデータ数が2の冪のものに限定されるが、任意の自然数Nに対しO(Nl ogN)の計算量でより高速に計算できる改良型アルゴリズムが存在する。

基本原理

1の原始N乗根のひとつをW_N(=exp(-2πi/N))とする。（W_Nを回転子と呼ぶ。）
データ数Nの列a=[f(0), f(1), f(2), ... ,f(N-2),f(N-1)]に対する通常のフーリエ変換は以下の式で表される。

F(u)=(1/N)Σ^N-1_n=0f(n)W_N^nu=F(f)(u)

つまり、離散フーリエ変換Fにより、数列a=[f(0), ... ,f(N-1)]が数列A=[F(0), ... ,F(N-1) ]に移されるということである。

F : a　→　A

この計算では、N個のデータaからN個のデータAを作るために、(N回の積+N回の和)×N回≒O(N²)の計算が要求される。

高速フーリエ変換では半分の長さの配列に分割していく作業をする。1回半分にすると(N/2回の積+N/2回の和)×N/2回×2であり一度分割すると計算回数はおよそ半分になる。最終的に1配列あたりO(log₂N)の計算量と、分割する作業がN-1回なので、全体としてO(Nl og₂N)程度になる。

高速化の準備

ここで、uを偶数と奇数の場合に場合分けしてみる。2N'=Nと置く。

u=2u'の時、

N×F(2u')=Σ^N-1_n=0f(n)W_N^n×2u'
　=Σ^N'-1_n=0f(n)W_2N'^2nu'+Σ^2N'-1_n=N'f(n)W_2N'^2nu'　=Σ^N'-1_n=0f(n)W_N'^nu'+Σ^N'-1_n=0f(n+N')W_N'^(n+N')u'　=Σ^N'-1_n=0f(n)W_N'^nu'+Σ^N'-1_n=0f(n+N')W_N'^N'u'W_N'^nu'　=Σ^N'-1_n=0(f(n)+f(n+N'))W_N'^nu'

u=2u'+1の時、

N×F(2u'+1)=Σ^N-1_n=0f(n)W_N^n(2u'+1)
　=Σ^N'-1_n=0f(n)W_2N'^n(2u'+1)+Σ^2N'-1_n=N'f(n)W_2N'^n(2u'+1)　=Σ^N'-1_n=0f(n)W_2N'^2nu'W_2N'²ⁿ+Σ^N'-1_n=0f(n+N')W_2N'^{(n+N')(2u'+1)}　=Σ^N'-1_n=0f(n)W_N'^nu'W_Nⁿ+Σ^N'-1_n=0f(n+N')W_2N'^2nu'W_2N'^2N'u'W_2N'ⁿW_2N'^N'　=Σ^N'-1_n=0f(n)W_N'ⁿW_N'^nu'+Σ^N'-1_n=0f(n+N')W_N'^nu'×1×W_N'ⁿ×(-1)
　=Σ^N'-1_n=0{(f(n)-f(n+N'))W_N'ⁿ}W_N'^nu

ここで、

f'₀(n)=f(n)+f(n+N')
f'₁(n)=(f(n)-f(n+N'))W_N'ⁿ

とおくと、

F(2u')=(1/2N')Σ^N'-1_n=0f'₀(n)W_N'^nu'= (1/2)F(f'₀)(u')
F(2u'+1)=(1/2N')Σ^N'-1_n=0f'₁(n)W_N'^nu'= (1/2)F(f'₁)(u')

となり、　奇数、偶数の場合それぞれにFを作用させることができる。
つまり、元の長さNの配列aを長さN/2=N'のa'₀とa'₁に合成＆分解し、その後それぞれにフーリエ変換を施せばよい。Nが2で割り切れるかぎり、この作業を再帰的に適用させる事ができる。

Nが2の冪であるときは最終的にフーリエ変換により、F:[f(0), f(1)]→[(f(0)+f(1))/2, (f(0)-f(1))/2]という計算を生成された配列の数だけ繰り返すこととなる。

アルゴリズムまとめ

M>0、N>1の自然数とする。

長さN=2^MLの配列aから長さN'=2^M-1Lの配列a₀, a₁を生成する作用素をS_M-1と置く。S_M-1:a→a₀, a₁。　各成分は、

f₀(n)=f(n)+f(n+N')
f₁(n)=(f(n)-f(n+N'))W_Nⁿである。

aにフーリエ変換F_M-1を施した配列Aの各成分は、偶数、奇数で分類することができ、2N'=Nとすると、
F(2u')=(1/2N')Σ^N'-1_n=0(f(n)+f(n+N'))W_N'^nu'=Σ^N'-1_n=0f₀(n)W_N'^nu'=(1/2)F(f₀)(u')
F(2u'+1)=(1/2N')Σ^N-1_n=0((f(n)-f(n+N'))W_Nⁿ)W_N'^nu'=Σ^N-1_n=0f₁(n)W_N'^nu'=(1/2)F(f₁)(u')
となる。

長さ2^MNの配列Aを上記偶数成分、奇数成分で並び替える作用素をT_M-1と置く。T_M-1:A→A。各成分は、
F(u)=F(2u')
F(u+N’)=F(2u'+1)
とする。

アルゴリズム全体では以下のような再帰的なチャートになる。

→	aを受け取る　　→ ↓	長さが2（または途中で打ち切り）→	aにフーリエ変換
↑	aの長さが2で割りきれる ↓		↓
↑	a₀, a₁を生成 ↓		Aを得て終了
↑	←a₀、a₁でスタートに戻る

例として16個のデータからなる配列に対するアルゴリズムを図示する。

元デ｜タ	f_-(0) f_-(1) f_-(2) f_-(3) f_-(4) f_-(5) f_-(6) f_-(7) f_-(8) f_-(9) f_-(10) f_-(11) f_-(12) f_-(13) f_-(14) f_-(15)	→ S₂	f₀(0) f₀(1) f₀(2) f₀(3) f₀(4) f₀(5) f₀(6) f₀(7) f₁(0) f₁(1) f₁(2) f₁(3) f₁(4) f₁(5) f₁(6) f₁(7)	→ S₁	f_0,0(0) f_0,0(1) f_0,0(2) f_0,0(3) f_0,1(0) f_0,1(1) f_0,1(2) f_0,1(3) f_1,0(0) f_1,0(1) f_1,0(2) f_1,0(3) f_1,1(0) f_1,1(1) f_1,1(2) f_1,1(3)	→ S₀	f_0,0,0(0) f_0,0,0(1) f_0,0,1(0) f_0,0,1(1) f_0,1,0(0) f_0,1,0(1) f_0,1,1(0) f_0,1,1(1) f_1,0,0(0) f_1,0,0(1) f_1,0,1(0) f_1,0,1(1) f_1,1,0(0) f_1,1,0(1) f_1,1,1(0) f_1,1,1(1)
	↓F₃		↓F₂		↓F₁		↓F₀
変換後	F(0) F(1) F(2) F(3) F(4) F(5) F(6) F(7) F(8) F(9) F(10) F(11) F(12) F(13) F(14) F(15)	→ T₂	F(0) F(2) F(4) F(6) F(8) F(10) F(12) F(14) F(1) F(3) F(5) F(7) F(9) F(11) F(13) F(15)	→ T₁	F(0) F(4) F(8) F(12) F(2) F(6) F(10) F(14) F(1) F(5) F(9) F(13) F(3) F(7) F(11) F(15)	→ T₀	F(0) F(8) F(4) F(12) F(2) F(10) F(6) F(14) F(1) F(9) F(5) F(13) F(3) F(11) F(7) F(15)

f_{m0, m1, m2}(0)とf_{m0, m1, m2}(１)をあわせると、f(0)からf(15)までの成分をすべて1度ずつ（位相をずらしながら）含んでいる。この時の和と積の取り方を表したグラフが蝶の羽のように見えることから、バタフライ演算と呼ばれたりする。

2から一般の素数pに拡張

Nが2の冪だけでなく、pの冪であるとき、式は以下のようになる。N'p=Nと置く。0≦m<p-1。

F(pu'+m)=Σ^N'-1_n=0(Σ^p-1_k=0W_p^kW_N^mnf(n+kN'))W_N'^u'n

上記の式に沿ってaからa₀,a₁,…,a_p-1を生成し、データ数がpで割り切れなくなった時点で通常のフーリエ変換を施す、というアルゴリズムになる。
解説したのはp=2のパターン。p=4の場合が最も高速になる様子。

ビット反転

表からDFTとFFTでは得られる配列F(u)の順番が異なることがわかる。

aにSを作用させて生成した配列の添え字とかっこ内の数字を見てみると、たとえば（0番目から数えて）11番目の要素f_1,0,1(1)は数字が1,0,1,1と並んでおり、2進数で1011は10進数で11となる。6番目の要素f_0,1,1(0)は0110であり、これは10進数で6である。
つまり、n番目の要素は2進数でnと一致するように0,1が並んでいるのである。

一方、AにTを作用させて並び替えた配列は、11番目はF(13)、6番目はF(6)である。実は、13は2進数で1101であり、ちょうど11と逆転している。他のすべての数も、2番目の要素は2進数表記を逆順にした数字になっている。

以下のデータ数16の場合の表を見るとわかりやすい。

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

→
2
進
数
化

0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

→
ビ
ッ
ト
反
転

0000
1000
0100
1100
0010
1010
0110
1110
0001
1001
0101
1101
0011
1011
0111
1111

→
10
進
数
化

0
8
4
12
2
10
6
14
1
9
5
13
3
11
7
15

この変換をビット反転と呼ぶ。FFTで計算した後には元の並びと一致させるためAにビット反転を施してやる必要がある。

実装

FFTは複素フーリエ変換なので複素数を扱える言語で無ければ実数からなるデータは実数列1つとすべてが0の配列の2つを渡すことになる。

エクセル VBAの例：長さNの配列Re(), Im()の2つをFFT ()に渡し、フーリエ変換させて書き換える。

実際には安全性のためNが2の冪かどうかなどを確認する行をはさむ必要がある。

'=======================================メイン=========================
Function FFT(Re(), Im(), N)
Call FFTsub(Re(), Im(), N, 0)
Call BitInversion(Re(), N)
Call BitInversion(Im(), N)
For I = 0 to N - 1
Re(I) = Re(I) / N　　　'最後にまとめて係数をかける
Im(I) = Im(I) / N　　　
Next
End Function
'=======================================FFTアルゴリズム====================
Function FFTsub(Re(), Im(), N, St)  ' Stは分割された配列の先頭
Theta = 2 * 3.14159265358 / N
If N >1 Then
     For I=0 To N / 2 - 1
          A0 = Re(I + St) + Re(I + St + N / 2)
          A1 = Im(I + St) + Im(I + St + N / 2)
          B0 = Re(I + St) - Re(I + St + N / 2)
          B1 = Im(I + St) - Im(I + St + N / 2)
          Re(I + St) = A0
          Im(I + St) = A1
          Re(I + St + N / 2) = B0 * Cos(Theta * I) + B1 * Sin(Theta * I)
          Im(I + St + N / 2) = B1 * Cos(Theta * I) - B0 * Sin(Theta * I)
     Next I
     Call FFTsub(Re(),Im(),N / 2, St)           '再帰呼び出し
     Call FFTsub(Re(),Im(),N / 2 , St + N / 2)       '再帰呼び出し
End If
End Function
'=======================================配列のビット反転======================
Function BitInversion(F(), N)
Dim Ary() As Variant
ReDim Ary(N)
For I = 0 To N - 1
     Ary(I) = F(I)
Next I
For I = 0 To N - 1
     F(I) = Ary(BIsub(I, N))
Next I
End Function
'========================================N以下の自然数をビット反転=============
Function BIsub(Num, N)
A = N - 1
B = 0
C = Num
While A > 0
     A = A \ 2
     B = B * 2 + (C Mod 2)
     C = C \ 2
Wend
BIsub = B
End Function

逆FFT

フーリエ変換と同様にして、逆フーリエ変換も定義できる。逆変換は、
f(n)=Σ^N-1_u=0F(u)W_-N^un=F'(F)(n)
と定義されるので、

f(2n')=Σ^N-1_u=0(F(u)+F(u+N'))W_-N'^un'=F'(F₀)(n)
f(2n'+1)=Σ^N-1_u=0((F(u)-F(u+N'))W_-N'^u)W_-N'^un'=F'(F₁)(n)

となる。計算の流れもほぼ同じ。したがって、プログラムは一部だけ変形させ、ほぼそのまま流用できる。

以下のプログラムでは赤色の部分を書き換えている。

'=======================================メイン=========================
Function FFTI(Re(), Im(), N)       '名前変更
Call FFTIsub(Re(), Im(), N, 0)　　　
Call BitInversion(Re(), N)
Call BitInversion(Im(), N)
　　　　　　　　　　　　　　　　　　'最後に掛ける係数を1にする
End Function
'===================================逆FFTアルゴリズム====================
Function FFTIsub(Re(), Im(), N, St) 　'名前変更 
Theta = 2 * 3.14159265358 / N
If N >1 Then
     For I=0 To N / 2 - 1
          A0 = Re(I + St) + Re(I + St + N / 2)
          A1 = Im(I + St) + Im(I + St + N / 2)
          B0 = Re(I + St) - Re(I + St + N / 2)
          B1 = Im(I + St) - Im(I + St + N / 2)
          Re(I + St) = A0
          Im(I + St) = A1
          Re(I + St + N / 2) = B0 * Cos(Theta * I) - B1 * Sin(Theta * I) '±逆転 
          Im(I + St + N / 2) = B1 * Cos(Theta * I) + B0 * Sin(Theta * I) '±逆転
     Next I
     Call FFTIsub(Re(),Im(),N / 2, St)         　　'呼び出し関数変更
     Call FFTIsub(Re(),Im(),N / 2 , St + N / 2)    '呼び出し関数変更   
End If
End Function

高速フーリエ変換単語

概要