◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇ 
 
高速フーリエ変換（FFT）の原理 
 
発行：エスオーエル株式会社 
https://www.sol-j.co.jp/ 
 
連載「X線CTで高精度寸法測定！？」 
2014年12月10日号　VOL.043 
 
平素は格別のお引き立てを賜り、厚く御礼申し上げます。 
X線CTスキャンによる精密測定やアプリケーション開発情報などをテーマに、 
無料にてメールマガジンを配信いたしております。 
 
◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇ 

 
 
最近、息子がよく「テューター、テューター」と言います。 
どうやら「機関車トーマス」のことのようです。 
トーマスの仲間の「パーシー」や「トビー」はちゃんと言えるのに、 
肝心の「トーマス」はまだ発音できないようです。 
 
 
さて、測定機の重要な基礎の一つにフーリエ変換があることは、 
何度もお伝えしてきました。 
 
X線CTにしても、干渉計のフリンジスキャンにしても、 
フーリエ変換なしでは、原理に立ち返って開発や提案をすることができません。 
装置設置場所の振動などの外乱を調査するにも、 
測定で得られた生データから必要な情報を取り出すフィルタリングを実行するにしても 
フーリエ変換が必須です。 
 
 
フーリエ変換を実際にコンピュータに計算させるには、 
DFT（離散フーリエ変換）を実行することになります。 
 
DFTは、W = exp(-2πi/N) と書いておくと、 
 
　　F(k) = Σ^{N-1}_{j=0} f(j) W^(kj) 
 
を計算せよということになります。 
 
この計算をそのまま実行すると、 
j=0～N まで f(j) と W^(kj) を N回掛け算することになります。 
N回掛け算することで、ようやくフーリエ変換後の1点が決まります。 
 
k=0～N まで全ての F(k) を求めるには、N×N 回の掛け算が必要です。 
 
足し算の回数も数えておくと、N(N-1) 回あります。 
 
 
では、高速フーリエ変換（FFT）と呼ばれる方法で、 
どれだけ計算量が減るかを見るとともに、 
FFTの基本原理をご紹介致します。 
 
まず、F(k) について、和を偶数番目と奇数番目の項に分けます。 
 
F(k) = Σ^{N/2-1}_{j=0} f(2j) W^(2jk) 
　　　　+ Σ^{N/2-1}_{j=0} f(2j+1) W^((2j+1)k) 
 
　　 = Σ^{N/2-1}_{j=0} f(2j) exp(-2πikj/(N/2)) 
　　　　+ Σ^{N/2-1}_{j=0} f(2j+1) exp(-2πik/N) exp(-2πikj/(N/2)) 
 
ここで、 
 
　　F[e](k) = Σ^{N/2-1}_{j=0} f(2j) exp(-2πikj/(N/2)) , 
　　F[o](k) = Σ^{N/2-1}_{j=0} f(2j+1) exp(-2πikj/(N/2)) , 
　　W[N] = exp(-2πi/N) 
 
と書くことにすると、 
 
F(k) = F[e](k) + W[N]^k F[o](k) 
 
と分割できることが分かります。この分割をさらに続けます。 
F[e](k) と F[o](k) を同様に分割することで、 
 
F(k) = F[ee](k) + W[N/2]^k F[eo](k) + W[N]^k ( F[oe](k) + W[N/2]^k F[oo](k) ) 
 
とすることができます。中身を書いておくと、 
 
　　F[ee](k) = Σ^{N/4-1}_{j=0} f(4j) exp(-2πikj/(N/4)) , 
　　F[eo](k) = Σ^{N/4-1}_{j=0} f(4j+2) exp(-2πikj/(N/4)) , 
　　F[oe](k) = Σ^{N/4-1}_{j=0} f(4j+1) exp(-2πikj/(N/4)) , 
　　F[oo](k) = Σ^{N/4-1}_{j=0} f(4j+3) exp(-2πikj/(N/4)) , 
　　W[N/2] = exp(-2πi/(N/2)) 
 
です。この調子で続けていくと、Σ^{0}_{j=0} にたどり着きます。 
そうすると、f(jの式) に j=0 を代入して、 
さらに exp(0) の部分は 1 になりますので、単純な式になります。 
 
例えば、N=2 のときは、 
　　F(k) = f(0) + W[2]^k f(1) 
 
N= 2^2 = 4 のときは、 
　　F(k) = f(0) + W[4]^k f(2) + W[2]^k f(1) + W[2]^k W[4]^k f(3) 
 
N= 2^3 = 8 のときは、 
　　F(k) = f(0) + W[8]^k f(4) + W[4]^k f(2) + W[4]^k W[8]^k f(6) 
　　　　　　+ W[2]^k f(1) + W[2]^k W[8]^k f(5) + W[2]^k W[4]^k f(3) 
　　　　　　+ W[2]^k W[4]^k W[8]^k f(7) 
 
と書けます。 
 
一般に、N = 2^m と書くと、F(k) は m回（つまり、log_2 N 回）分割できます。 
そして、1回の分割で W[]^k という因子が全体の半分に掛かりますので、 
m回分割後には、1つの項当たりに m/2 回の掛け算が実行されてる換算になります。 
 
その項が N個できますので、掛け算の回数は、N×m/2 , 
つまり、(N/2) log_2 N となります。 
 
ちなみに足し算は、1つの k につき、分割した回数（m = log_2 N 回）必要なので、 
k=0～N の全てで、N×log_2 N 回の足し算ということになります。 
 
 
計算量が分かったところで、どれ位早くなるのかを比較します。 
 
N×N 回のオーダーと N×log_2 N 回のオーダーは、どの位違うのでしょうか。 
N=1024×1024 の計算をさせたとします。 
 
その場合、 
　　N×N = 1024^4 = 2^40 = 1099511627776 , 
　　N×log_2 N = 1024^2 × 20 = 20971520 
を比較し、52428.8 倍の違いがあることが分かります。 
 
もし、FFTを使うと1秒で完了する計算だったとすると、 
FFTなしでは、14時間半も待たされることになります。 
 
 
FFTの基本原理は、上に述べた通りですが、 
実用的なアルゴリズムの実装には、もうひと手間必要です。 
 
機会があれば、バタフライ演算やビット逆転についてご紹介したいと思います。 
 
また、FFTの分割を最後まで実行せずに、N=4 までで止めておいて、 
sin や cos の値が 0, ±1 となる性質を利用して、 
掛け算をさらに省略する方法などもあり、FFTの奥は深いです。 


--
高野智暢


☆TomoScope専門サイトはこちら☆