高速フーリエ変換単語

コウソクフーリエヘンカン
  • 3
  • 0pt
掲示板へ

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

もしかしてファイナルファンタジータクティクスFFT

概要

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

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

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

基本原理

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

F(u)=ΣN-1n=0f(n)WNnu=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(N2)の計算が要される。

高速化の準備

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

u=2u'の時、

F(2u')=ΣN-1n=0f(n)WNn×2u'
 =ΣN'-1n=0f(n)WN2nu'+Σ2N'-1n=N'f(n)WN2nu'
 =ΣN'-1n=0f(n)WN'nu'+ΣN'-1n=0f(n+N')WN'(n+N')u'
 =ΣN'-1n=0f(n)WN'nu'+ΣN'-1n=0f(n+N')WN'N'u'WN'nu'
 =ΣN'-1n=0(f(n)+f(n+N'))WN'nu'

u=2u'+1の時、

F(2u'+1)=ΣN-1n=0f(n)WNn(2u'+1)
 =ΣN'-1n=0f(n)WNn(2u'+1)+Σ2N'-1n=N'f(n)WNn(2u'+1)
 =ΣN'-1n=0f(n)WNnWN2nu'+ΣN'-1n=0f(n+N')WN(n+N')(2u'+1)
 =ΣN'-1n=0f(n)WNnWN'nu'+ΣN'-1n=0f(n+N')WN2nu'WN2N'u'WNnWNN'
 =ΣN'-1n=0f(n)WNnWN'nu'+ΣN'-1n=0f(n+N')WN2nu'×1×WNn×(-1)
 =ΣN'-1n=0{(f(n)-f(n+N'))WNn}WN'nu

ここで、

f'0(n)=f(n)+f(n+N')
f'1(n)=(f(n)-f(n+N'))WNn

とおくと、

F(2u')=ΣN'-1n=0f'0(n)WN'nu'= F(f'0)(u')
F(2u'+1)=ΣN'-1n=0f'1(n)WN'nu'= F(f'1)(u')

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

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

アルゴリズムまとめ

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

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

f0(n)=f(n)+f(n+N')}
f1(n)=(f(n)-f(n+N'))WNn
である。

aにフーリエ変換FM-1を施した配列Aの各成分は、偶数奇数で分類することができ、2N'=Nとすると、
F(2u')=ΣN'-1n=0(f(n)+f(n+N'))WN'nu'=ΣN'-1n=0f0(n)WN'nu'=F(f0)(u')
F(2u'+1)=ΣN-1n=0((f(n)-f(n+N'))WNn)WN'nu'=ΣN-1n=0f1(n)WN'nu'=F(f1)(u')
となる。

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

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

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

例として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)









S2
f0(0)
f0(1)
f0(2)
f0(3)
f0(4)
f0(5)
f0(6)
f0
(7)
f1(0)
f1(1)

f1(2)
f1(3)
f1(4)
f1(5)
f1(6)
f1(7)









S1
f0,0(0)
f0,0(1)
f0,0(2)
f0,0(3)
f0,1(0)
f0,1(1)
f0,1(2)
f0,1
(3)
f1,0(0)
f1,0(1)

f1,0(2)
f1,0(3)
f1,1(0)
f1,1(1)
f1,1(2)
f1,1(3)









S0
f0,0,0(0)
f0,0,0(1)
f0,0,1(0)
f0,0,1(1)
f0,1,0(0)
f0,1,0(1)
f0,1,1(0)
f0,1,1
(1)
f1,0,0(0)
f1,0,0(1)

f1,0,1(0)
f1,0,1(1)
f1,1,0(0)
f1,1,0(1)
f1,1,1(0)
f1,1,1(1)
F3 F2 F1 F0









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)








T2
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)








T1
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)








T0
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)

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

2からpに拡張

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

F(pu'+m)=ΣN'-1n=0(Σp-1k=0WpkWNmnf(n+kN'))WN'u'n

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

ビット反転

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

aにSを作用させて生成した配列の添え字とかっこ内の数字を見てみると、たとえば(0番から数えて)11番の要素f1,0,1(1)は数字が1,0,1,1と並んでおり、2進数で1011は10進数で11となる。6番の要素f0,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)
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)=(1/N)ΣN-1u=0F(u)W-Nun=F'(F)(n)
と定義されるので、

f(2n')=ΣN-1u=0(F(u)+F(u+N'))W-N'un'=F'(F0)(n)
f(2n'+1)=ΣN-1u=0((F(u)-F(u+N'))W-Nu)W-Nun'=F'(F1)(n)

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

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

'=======================================メイン=========================
Function FFTI(Re(), Im(), N)
Call FFTIsub(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 I
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

関連動画

関連静画

関連商品

関連項目

【スポンサーリンク】

  • 3
  • 0pt
スマホ版URL:
https://dic.nicovideo.jp/t/a/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B

この記事の掲示板に最近描かれたお絵カキコ

お絵カキコがありません

この記事の掲示板に最近投稿されたピコカキコ

ピコカキコがありません

高速フーリエ変換

まだ掲示板に書き込みがありません…以下のようなことを書き込んでもらえると嬉しいでーす!

  • 記事を編集した人の応援(応援されると喜びます)
  • 記事に追加して欲しい動画・商品・記述についての情報提供(具体的だと嬉しいです)
  • 高速フーリエ変換についての雑談(ダラダラとゆるい感じで)

書き込みを行うには、niconicoのアカウントが必要です!