連載「X線CTで高精度寸法測定!?」

◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇ 高速フーリエ変換(FFT)の原理 発行:エスオーエル株式会社 http://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専門サイトはこちら☆ ●┳┳┳●━━━━ 連 絡 先 ━━━━━━━━━━━━━ ┣╋╋○ エスオーエル株式会社 ( SOL ) ┣╋○ 〒335-0012 埼玉県戸田市中町1-34-1 ┣○ Tel: 048-441-1133 Fax: 048-445-1678 ● Email: sales@sol-j.co.jp Web: http://www.sol-j.co.jp    --デモ測定を承ります-- 詳細は上記Webサイトまで

測定機に関する理論や解析技術の情報発信として、定期的にメールマガジンを無料配信しています。メルマガにご登録頂くと、最新情報やWEB未公開の情報をいち早くお手元にお届けすることができます。

メルマガご登録はこちらから