図1: 時間間引き型 FFT の概念図
by 池田 幹男 since 2002 Oct 22
このページの内容は、 四日市大学教育支援システム のコース Cooley Tukey 型高速フーリエ変換に移動しました。内容はゲストでログインして見ることができます。
プログラムはこのアーカイブを参照して下さい。
FFT (Fast Fourier Transform)とは DFT(Discrete Fourier Transform) の高 速算法(アルゴリズム)のことで、高速フーリエ変換と呼ぶ。DFT は式(1) のように定義される。
F(k) = Σn=0N-1 s(n) Wnk/N (n,k=0,1,...,N-1) ....(1)
ただし、Wnk/N = exp(j2πnk/N) = cos(2πnk/N) + j sin(2πnk/N) であり、 j は虚数単位である。
これを素朴な方法で計算すると、N2 回の乗算と加算が必 要になる。N が小さい場合には、十分実用的な速度で計算可能である が、大きくなると非常に遅くなる。
ここで、N が MとKを約数に持つとすると (1)式は (2) 式に変換できる。
F(pK+k) = Σq=0K-1 Σm=0M-1 s(qM+m) W(pK+k)(qM+m)/N (2)
順次変形して、(5)式が導き出される。
F(pK+k) = Σq=0K-1 Σm=0M-1 s(qM+m) W(pqKM +pKm +qkM+km)/N (3)
F(pK+k) = Σq=0K-1 Σm=0M-1 s(qM+m) W(pKm+qkM+km)/N (4)
F(pK+k) = Σq=0K-1 Σm=0M-1 s(qM+m) Wpm/M Wqk/K Wkm/N (5)
(5)式を見てみると、Σq=0K-1 s(qM+m) Wqk/K の部分が、信号 s(n) を M 個おきにとったものの K 点の DFT となっている。この部分 DFT を Gm(k) (m=0,...,M-1) とすると、(6)式となり、
F(pK+k) = Σm=0M-1 Gm(k) Wkm/N Wpm/M (6)
F(pK+k) (p=0,1,...,M-1) は Gm(k) (m=0,1,...,M-1)に Wkm/N を乗じたものの M 点の DFT となっている。 この操作を図1 にしめす。図1をみると、s(n) を M の 剰余系で分類したものを、K点の DFT した後に、K個おきにとっ たものを M 点の DFT することになる。この操作の場合にはまず、信号を M の剰余系毎に整列しなおすことによって、周波数順序の揃った DFT の結果を得ることができる。これが、時間間引き FFTの考え方であ る。(この図や、(6)式は、二次元 DFT を行なっている。従って、ちょっとし たプログラムの変更で二次元 DFT も容易に実現可能である。多次元 DFT もマ ルチステージの FFT をもとに可能となりそうである。)
図1: 時間間引き型 FFT の概念図
同じ操作を時間順序が揃った順に並べ直したものを図2 に示す。図2をみると、 s(n) を M 個おきとったものを、K点の DFT し た後に、M点の DFT することになる。この操作の場合には出力のスペクトルを M の剰余系によって、整列しなおすことによって、周波数順序の揃っ た DFT の結果を得ることができる。これが、周波数間引き FFTの考 え方である。
図2 : 周波数間引き型 FFT の概念図
さらに、K,M が別の約数を持つならば、部分 DFT も同じような方法で約数の DFT に分解して計算することができる。次の節では FFT でどのように計算量が 削減されるかを見てみよう。
上記より、N = K1 K2 ... KP と 約数Kp (以下ではこの Kp のことを 基数とも呼ぶ)に分解できるとすると、必要な計算量は、各 ステージで Kp点の部分 DFT を N/Kp 回計算し、Wkm/N を乗ずる操作が N 回あり、最終 段では Wkm/N を乗ずる必要がないことから、 Kp点の DFT を素朴な方法で実装しても、
必要となると推測される。複素数の乗算は実数の乗算が 4 回、加算が 2 回必 要であるとするとして(乗算を減らすことはできるが、加算が増えるので、こ こでは考慮しない。)、さらに Wpm/M の対称性や 0, 1 と なる部分を利用し、K が素数の(約数を持たない)場合、
乗算が
4N(P-1) + Σp=1P
(N/Kp)(Kp-1)2
=
2NP-4N + Σp=1P
NKp
加算が
2N(P-1) + Σp=1P
(N/Kp)(Kp2 +
2Kp-1)
= 4NP-2N +
Σp=1P NKp
回程度となる。 Kp がすべて等しい K と仮定すると、 P=logKN となるので、計算量は、
となる。K=2 を代入すると、合計 10Nlog2N-6N となるが、実際の K=2 の場合は、もう少し効率化が可能で、乗算が、2N log2N - 4N 回程度、加算が、3N log2N - 2N 回程度になる。加算と乗算のウェイト が同じならば、合計 5Nlog2N - 6N となる。 同様に一般のKの場合、N(6+2K)log2N /log2K - 6N となる計算となる。従って、オーダーとしては、計算量は同じで あり、1/log2K 倍されるので、多少大きな約数となっても、 効率がさほど大きく落ちるわけではない。例えば、K=17としても、 9.79Nlog2N - 6N であり、基数 2 の場合 と比較しても 倍ほどは遅くならない。面白いのは、この式による推定値では、 K=5 のとき最小値、約6.89Nlog2N を持つ ことである。実際には、小さな K の場合には、乗算を減らすことがで きるので、この式の予測値よりも小さな値を持つ。従って、Nに比べて 十分小さな約数をもつ FFT の計算量は、最良のばあいで、 4Nlog2N (基数 8,16)程度であり、基数2の場合と 比較して 20 % 程度しか効率化できないことになる。通常は (5〜 10)Nlog2N 程度であると推測される。表 1 に基数 が小さく、Nが十分に大きい場合の FFT の計算量を示す。実際には、 Wnk/N が 1,j,-1,-j になるものが相当量 含まれていて、この部分の計算を省くことができるので、オーダーには影響し ないが、もう少し計算量は小さくなる。
実際には、素数の点数に対する高速アルゴリズムも存在するが、ここではまだ 考慮していない。
このように、素朴な DFT を利用しても、point数が小さな約数を持てば、 FFT のアルゴリズムを利用して十分に高速化が可能であることがわかる。しかも、 FFT を用いて計算した方が、素朴な DFT を使用した場合よりも精度は向上す る。従って、N がある程度よりも大きい場合には FFT を利用した方が圧倒的 に有利である。
| 乗算の回数 | 加算の回数 | 合計 | |
|---|---|---|---|
| K=2 | 2Nlog2N - 4N | 3Nlog2N - 2N | 5Nlog2N - 6N |
| K=3 | 3.36Nlog2N - 4N | 2.52Nlog2N - 2N | 5.88Nlog2N - 6N |
| K=4 | 1.50Nlog2N - 4N | 2.75Nlog2N - 2N | 4.25Nlog2N - 6N |
| K=5 | 3.44Nlog2N - 4N | 2.76Nlog2N - 2N | 6.20Nlog2N - 6N |
| K=6 | 3.22Nlog2N - 4N | 1.55Nlog2N - 2N | 4.77Nlog2N - 6N |
| K=7 | 3.66Nlog2N - 4N | 3.05Nlog2N - 2N | 6.72Nlog2N - 6N |
| K=8 | 1.33Nlog2N - 4N | 2.75Nlog2N - 2N | 4.08Nlog2N - 6N |
| K=17の推定値 | 4.65Nlog2N - 4N | 5.13Nlog2N - 2N | 9.79Nlog2N - 6N |
| 一般の K の推定値 | N(2+K)log2N/log2K - 4N | N(4+K)log2N/log2K - 2N | N(6+2K)log2N/log2K - 6N |
| 素朴な DFT の場合 | N2 | N2 | 2N 2 |
Kが大きな素数の場合、乗算は、K2-2K+1 回 加算は、K2+2K-1回なので、合計 2K2 回。
ちなみに、より素朴な方法で計算すると、乗算は 4(K-1)2 回 加算はほぼ 4K2 回の合計 8K2-12K+3 回となり、ちょっとした工夫で、計算量を1/4にで きる。これは、大きな素数は必ず奇数であるので、奇数-1は偶数である。従っ て、素数点数の FFTからわかるように、少なく とも一つの約数 2 によって高速化可能であるのと等価である。
N=KM とすると、FFT の算法では、 2MK2+2KM2+6N回の計算が必要 となる。2N(K+M)+6N であり、計算量はこれだけでほぼ (K+M)/N に削減できる。K=Mとすると、 (K+M)/N=sqrt(2)/sqrt(N) であり、これだけで N2 の計算量を N1.5にできる。ただし、 2しか約数をもたない場合には、(K+M)/N はぼぼ 1/2となり、 半分にしか減らすことができない。
上記の予測では基数を 8 にしても、基数 2の場合と比較して、さほど FFT は高速化されないことになる。しかし、実際のコンピュータでは基数を8とす ると、かなりの高速化が可能である。(場合にもよるが、倍程度は高速になる。) これは、基数 2 の場合は、基数 8 の場合と比較すると、3倍も同じメモリを 走査する必要となり、キャッシュの入れ替えが頻繁に生じるためであると考え られる。実際問題、メモリのアクセス順序を変更するだけでも高速化が可能で ある。さらに、ループのオーバーヘッドも基数2の場合の方が相対的に多くな り、遅い分岐命令を頻繁に実行することになる。したがって、実際には多少の 複雑さを犠牲にして、基数の大きな FFT のほうが有利になる可能性が高い(実 際問題、素数基数の FFT でさえ、基数2よりも高速化できる可能性さえある。 実際のプロファイルでも、基数 7 の FFT は意外と速い。) メモリが余分に必 要となるが、sin, cos を予め計算してテーブルとして持っていると、メモリ アクセス順序の選択方法が柔軟にできるようになるので、単に sin, cos の 計算を省くだけよりも高速化が可能である。
実際の速度の問題はさらに FFT の高速化を複雑にする。プロセッサの種類や キャッシュのサイズによって、速度が大幅な影響を受けることになる。従って、 基数の組合せの選択は実際の実行速度のプロファイリングによって変化するこ とになる。
FFT を高速化するためには、
精度についてはここでは、詳しくは説明しない。有名な Oppenheim, Shafer の Digital Signal Processing ではかなり丁寧に説明してある。ここでは、 単に直観的に説明するが、FFT では、加算を木構造を利用したのと同じ効果が ある。もちろん素朴な DFT でも、木構造を利用すれば、精度を向上させるこ とが可能であるが、高速化も同時に達成することが可能な FFT のほうが有利 となる。
浮動小数点の誤差は、主として加算で生じる。小さな数をいっぱい加算するこ とを考えれば、すぐにわかる。例えば、1.0E-6 を 1.0E6 回加算する例を考え 見ればわかると思う。ある時点で、1.0E-6 の合計は、1.0E-6 よりも十分に大 きくなるが、その時の誤差は、加算した結果の LSB 程度であるとすると、相 対誤差は加算した結果の大きい方に依存することになる。つまり、和が大きく なれば、大きくなるほど加算する毎に、LSB 分の誤差が増加することになる。 (実際にはランダムウォークとみなす必要があるが) つまり、1.0E6 回の加算 を行なえば、最悪 6 桁精度が落ちると考えれば良い。積が入る場合には、加 算での誤差がそのままスケーリングされると考えればよい。
誤差を減らすためには、加算の結果が最小になるペアを探し出して加算し、加 算したもとの二個データを削除し、加算した結果をデータとして付け加えると いう手順を繰り返し、データが一個になるまで繰り返せば良い。(類似した方 法でのシミュレーションで相当精度を向上させることができることがわかって いる。二進木と比較しても、10000点で10進1桁以上の精度向上が可能である。) しかし、この方法ではデータ依存のうえに、余りに手間がかかるので、平均的 に同じ程度大きさのデータの場合には、単純に二進木で加算すれば良い。ここ まで書けば、もうだいたいわかると思うが、FFT では各ステージで木を構成し て加算しているのと同じ効果があることは、そのアルゴリズムからわかる。精 度劣化は、木の深さ、つまり FFT の段数に依存することになる。
このような誤差の問題は、デジタルフィルタで生じる。フィルタを通した結果、 雑音を発生させるという皮肉な結果になることさえある。畳み込みの演算は、 FFT で可能であるので、FFT を使用した方が精度面でも速度面でも安全である ことがわかる。意外なことではあるが、世の中には 60000 次のフィルタリン グをするアプリケーションが当たり前のように存在するのである。
世の中には、こんなページよりももっと役に立つページがいっぱいありま す。そちらを参照して下さい。