皆さん、こんにちは。
今回がフーリエ変換における数値解析の最終回で、前回紹介した離散フーリエ変換フーリエ変換を数値解析で ～離散フーリエ変換～ - ちょぴん先生の数学部屋 (hatenablog.com)をより高速に計算する「高速フーリエ変換」というアルゴリズムを紹介します。
※高速フーリエ変換は、その英訳 Fast Fourier Transformの頭文字を取って、しばしばFFTと略されて呼ばれます。
1.高速フーリエ変換
1-1. コンセプト
まず、離散フーリエ変換・逆変換の式を再掲します。
f_N(n) (n=0,1,2,・・・,N-1)はN個のデータからなるデジタル信号、F_N(k) (k=0,1,2,・・,N-1)はそのデジタル信号のk番目の離散スペクトルです。
そして、離散フーリエ変換を実行するDFT行列の定義も再掲します。
α_Nは絶対値1で偏角が-2π/Nの複素数なので、N乗すると1になり、α_Nの累乗はN個おきの周期性を持つことが分かります。
高速フーリエ変換は、α_Nのこの性質を利用してW_Nのサイズを半分ずつ減らしていく、というコンセプトのアルゴリズムになります。
1-2. アルゴリズムの中身
では、具体的なアルゴリズムについて説明していきます。
ここでは例としてN=4の場合について取り上げます。
この時、離散フーリエ変換は、成分表示すると
と書けます。
次の工程がテクニカルなのですが、⑤を「偶数行目だけでまとめたもの」と「奇数行目だけでまとめたもの」の2つに分けます。
具体的には次のような感じです。
まずは、偶数行目だけ集めた⑥に注目します。
であることに注意して、α_4の指数を下げてあげると、
係数行列の左側と右側が全く同じ行列になります。
よって、この2×4の係数行列は、実質2×2行列にすることが可能です。
ここで、この2×2行列の正体を見ておきます。
一般に
（訳：αの番号と指数を同時に半分にしても、値は変化しない）
が成り立つので、
係数行列は、
W_2そのものになります。
よって、⑨に注目した時、仮に
という読み替えをしたとすると、
⑨は、実質N=2の場合の離散フーリエ変換に帰着します。
同様に、奇数行目だけ集めた⑦についても見ていきます。
こちらは⑧に加えて
となることも利用して指数を下げると、
係数行列の右側は左側を-1倍しただけで実質同じ行列になっています。
よって、こちらも2×4行列を実質2×2行列にすることができます。
ここで、偶数行目の場合と同様に⑮の係数行列がW_2になってくれれば好都合ですね。そうなるように後ろのベクトルを細工すると次のようにできます。
よって、仮に
という読み替えをしたとすれば、
やはり実質N=2の場合の離散フーリエ変換に帰着できます。
以上の工程でN=4からN=2に半減できたので、最後のN=2からN=1への半減を続けていきます。
偶数行目からできた⑬と奇数行目からできた⑱は本質的に同じ式なので、説明の重複を避けるため⑬に絞って考えます。
⑬も偶数行目と奇数行目に分けます。
ここまでくればもう直接計算が可能なのでやってしまうと、
という感じになります。
このように、データN個の離散フーリエ変換を、データ数が半分のN/2個の離散フーリエ変換に変えていくのが、高速フーリエ変換の肝の部分であり、この一連の操作をサイズが1になるまで続けます。
このようなアルゴリズムの性質上、高速フーリエ変換が使えるのは、データ数Nが2進数のときだけです。むしろ、この事を見越してデータ数Nが2進数になるようにアナログ信号からサンプリングするのが一般的です。
一般化すると次のようになります。
1-3. 高速フーリエ変換の計算量
名前に「高速」と付いてるだけあって、ちゃんと通常の離散フーリエ変換に比べて計算量が減っています。
一般に、N個のデータからなるデジタル信号を（通常の方法で）離散フーリエ変換するには、前回触れた通りN^2回の計算が必要です。
高速フーリエ変換では、データの数を次々に半分にしていき、最終的には「1個のデータからなるデジタル信号を離散フーリエ変換する」をN回分行う計算に置き換えるので、計算量を見積もると、
N回となります。
そして、サイズを半減にする工程が、N=2^mの場合m回発生します。
よって、トータルの計算量は両者をかけたもの、つまり
N*logN回ということになります。元々の離散フーリエ変換がN^2回かかっていたのと比べると、確かに高速化できています。
その時短効果は前回触れた通りです。
2.例題
では、高速フーリエ変換を実際に例題を通してやってみましょう。
前回同様、次のようなN=4のデジタル信号を高速フーリエ変換してみます。
離散フーリエ変換の式を偶数行目と奇数行目に分けて、N=4をN=2に減らします。
この状態からさらにN=1に減らすと、(B),(C)それぞれで偶数行目、奇数行目を考えると、
となって、全ての離散スペクトルの値が求まります。
あとは、これをkの番号順に並べてあげれば、
前回と全く同じ結果が得られます。
高速フーリエ変換をPythonで実装すると次のようになります。
今回の例でプログラムを実行すると、
手計算の結果と一致します。
