Bluesteinのアルゴリズムとは (ブルースタインノアルゴリズムとは) [単語記事]

Bluesteinのアルゴリズムとは、任意の長さの離散データを高速フーリエ変換するためのFFT アルゴリズムである。

概要

ここでは、元のデータaを周期Nの周期的データと仮定する。

通常のフーリエ変換は任意の長さの配列を計算できるがO(N²)の計算量が必要であり、FFTはO(Nl ogN)で計算できるがNが2の冪に限定されるという問題がある。

Bluesteinのアルゴリズムは任意の長さの配列をO(Nl ogN)で計算できるアルゴリズムである。
原理は、フーリエ変換を畳み込みに変換して、畳み込みをFFTで計算するというものである。畳み込みは定義通りだとO(N²)の計算量になる。

配列の長さを2倍以上に拡張し、アルゴリズムの一部に2の冪のFFTを3回含むので、FFTと比べて6~8倍程度の計算時間になる。最適化しないうちは10倍以上はかかると想定しておいた方がよいだろう。

解説

W_N(m)=exp(πim/N)とする。
長さNのデータa=(a(0), a(1), ...,a(N-1))の離散フーリエ変換の定義は以下の通り。

F(u)=F(a)(u)=(1/N)Σ^N-1_n=0a(n)W_N(-2un)

(u-n)²=u²-2un+n²　であることに注意すると、2un=u²+n²-(u-n)²　なので、

W_N(-2un)＝W_N(-u²-n²+(u-n)²)
F(u)=(1/N)Σ^N-1_n=0a(n)W_N(-u²)W_N(-n²)W_N((u-n)²)
=(1/N)W_N(-u²)Σ^N-1_n=0a(n)W_N(-n²)W_N((u-n)²)

ここで、

a(m)W_N(-m²)=p(m)
W_N(m²)=q(m)

と置くと、

F(u)=q(-u)(1/N)Σ^N-1_n-0p(n)q(u-n) = q(-u)(p*q)(u) = q(-u)r(u)

となるので、a(n)のフーリエ変換をp(n)とq(n)の畳み込みに置き換えることができる。
p(n)とq(n)の畳み込みはフーリエ変換→各成分の積→フーリエ逆変換で置き換えることができる。
また、p(n)とq(n)を適切な構造の長さが2の冪の配列に拡張すれば、高速フーリエ変換を適用できるようになる。
q(n+N)=W_N(n²+N(2n+N))=q(n)なので、拡張されたp(n),q(n)は以下の通り。

p’(n) = p(n) = a(n)W_N(-n²) (0≦n<N)
p’(n) = 0　　　　　　(N≦n<2^M)
q’(n) = W_N(n²) 　　 (0≦n<2^M)

Mは、2^M-1<2N≦2^M となることが条件。

拡張後に得られた配列F'(u)は元の配列F(u)を含むので、F'(u)から適切な位置で長さNの分だけ切り取る。配列を長くすることが原因で畳み込みの振幅がN/M倍に変わるので最後にM/N倍する必要がある。

まとめ

全体のアルゴリズムは以下のようになる。k,mは自然数。

aを受け取る ↓	→	F’(u)=r(u)q(-u)を得る
aから長さ2^Mのp’(k),q’(k)を生成 ↓	↑	F'(u)からF(u)を取り出す ↓
p(k), q(k)をP(m), Q(m)に高速フーリエ変換 ↓	↑	振幅を整える ↓
R(m)=P(m)Q(m)を逆高速フーリエ変換	→↑	Aを得て終了

実装

VBAでは以下のようになる。データ部分の実数列Re()、全て0の構成された虚数列Im()を受け取り、直接書き換える。
2の冪のFFTと逆FFTは実装済みとする。

Function Bulestein(Re(), Im(), N)
Pi=3.1415926535
M = 1
While M < N
M = 2 * M
Wend
'==========================================2の冪には直接FFTを適用
If N = M Then      　　　
Call FFT(Re(), Im(), N)
'==========================================アルゴリズム開始
Else                     
M = 2 * M
Dim C() As Variant
ReDim C(M)
Dim S() As Variant
ReDim S(M)
Dim WC() As Variant
ReDim WC(M)
Dim WS() As Variant
ReDim WS(M)

For I = 0 To M - 1      '================================q(m)
WC(I) = Cos(Pi * I * I / N)
WS(I) = Sin(Pi * I * I / N)
Next I
For I = 0 To N - 1      '=================================p(m)
C(I) = Re(I) * WC(I) + Im(I) * WS(I)
S(I) = Im(I) * WC(I) - Re(I) * WS(I)
Next I
For I = N To M - 1
C(I) = 0
S(I) = 0
Next I
'================================================================畳み込み開始
Call FFT(C(), S(), M)
Call FFT(WC(), WS(), M)

For I - 0 to M - 1
X = C(I) * WC(I) - S(I) * WS(I)
Y = S(I) * WC(I) + C(I) * WS(I)
C(I) = X
S(I) = Y
Next I

Call FFTI(C(), S(), M)
'=================================================================畳み込み終了
Y = M / N  　'===========================================フーリエ変換に合わせる係数
For I = 0 To N - 1
X = I + N       '===========================================切り取る位置
Re(I) = (C(X) * Cos(Pi * X * X / N) + S(X) * Sin(Pi * X * X / N)) * Y
Im(I) = (S(X) * Cos(Pi * X * X / N) - C(X) * Sin(Pi * X * X / N)) * Y
Next I
End If
End Function

逆Bluesteinアルゴリズム

当然逆変換も存在する。フーリエ逆変換の定義から、

F^-1(u)=Σ^N-1_n=0a(n)W_N(2nu)=q(u)Σ^N-1_n=0p(n)q(u+n)

となり、畳み込みではなく相関関数を計算することになる。
全体としては係数と一か所だけ逆FFTに変えるだけで逆変換になる。

赤色が変更した部分。

Function BulesteinI(Re(), Im(), N)
Pi=3.1415926535
M = 1
While M < N
M = 2 * M
Wend
'==========================================2の冪には直接FFTを適用
If N = M Then      　　　
Call FFTI(Re(), Im(), N)　　　'===========逆FFT
'==========================================アルゴリズム開始
Else                     
M = 2 * M
Dim C() As Variant
ReDim C(M)
Dim S() As Variant
ReDim S(M)
Dim WC() As Variant
ReDim WC(M)
Dim WS() As Variant
ReDim WS(M)

For I = 0 To M - 1      '================================q(m)
WC(I) = Cos(Pi * I * I / N)
WS(I) = Sin(Pi * I * I / N)
Next I
For I = 0 To N - 1      '=================================p(m)
C(I) = Re(I) * WC(I) + Im(I) * WS(I)
S(I) = Im(I) * WC(I) - Re(I) * WS(I)
Next I
For I = N To M - 1
C(I) = 0
S(I) = 0
Next I
'================================================================畳み込み開始
Call FFTI(C(), S(), M)    '================================逆FFT
Call FFT(WC(), WS(), M)

For I - 0 to M - 1
X = C(I) * WC(I) - S(I) * WS(I)
Y = S(I) * WC(I) + C(I) * WS(I)
C(I) = X
S(I) = Y
Next I

Call FFTI(C(), S(), M)
'=================================================================畳み込み終了
Y = 1  　'===========================================フーリエ変換に合わせる係数
For I = 0 To N - 1
X = I      '===========================================切り取る位置
Re(I) = (C(X) * Cos(Pi * X * X / N) + S(X) * Sin(Pi * X * X / N)) * Y
Im(I) = (S(X) * Cos(Pi * X * X / N) - C(X) * Sin(Pi * X * X / N)) * Y
Next I
End If
End Function