高速フーリエ変換

高速フーリエ変換は、離散フーリエ変換を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換と呼ぶ。

1. 概要
-----

複素関数 //f//}} の離散フーリエ変換である複素関数 //F//}} は以下で定義される。

このとき、1=mset を標本点と言う。

これを直接計算したときの時間計算量は、ランダウの記号を用いて表現すると //O//}} である。

高速フーリエ変換は、この結果を、次数N}}が2の累乗のときに //O//}} の計算量で得るアルゴリズムである。より一般的には、次数が 1=//N// = ∏ //nsub と素因数分解できるとき、//O//(//N//∑//n//sub の計算量となる。次数が 2}} の累乗のときが最も高速に計算でき、アルゴリズムも単純になるので、0}} 詰めで次数を調整することもある。//

高速フーリエ変換を使って、畳み込み積分などの計算を高速に求めることができる。これも計算量を //O//}} から //O//}} まで落とせる。

現在は、初期の手法 をより高速化したアルゴリズムが使用されている。

1.1. 逆変換
--------

逆変換は正変換と同じと考えて良いが、指数の符号が逆であり、係数 1///N//}} が掛かる。

高速フーリエ変換のプログラム中、どの符号が逆転するかを一々分岐させると、分岐の判定に時間がかかり、パフォーマンスが落ちる。一方、正変換のプログラムと、逆変換のプログラムを両方用意しておくことも考えられるが、共通部分が多いため、無駄が多くなる。このため、複素共役を使った次のような方法が考えられる。

離散フーリエ変換を

で定義したとき、逆変換は

となる。

このため、//F//}} の離散フーリエ逆変換を求めるには、

 複素共役を取り、overline を求める、 overline の正変換の離散フーリエ変換を高速フーリエ変換で行う、 その結果の複素共役を取り、N}} で割るとすれば良く、正変換の高速フーリエ変換のプログラムがあれば、逆変換は容易に作ることができる。

2. アルゴリズム
---------

2.1. クーリー–テューキー型FFTアルゴリズム
-------------------------

クーリー–テューキー型アルゴリズムは、代表的な高速フーリエ変換  アルゴリズムである。

分割統治法を使ったアルゴリズムで、1=//N// = //N//sub のサイズの変換を、より小さいサイズである //N//sub, //N//sub のサイズの変換に分割していくことで高速化を図っている。

最もよく知られたクーリー–テューキー型アルゴリズムは、ステップごとに変換のサイズをサイズ //N///2}} の2つの変換に分割するので、2}} の累乗次数に限定される。しかし、一般的には次数は 2}} の累乗にはならないので、素因数が偶数と奇数とで別々のアルゴリズムに分岐する。

伝統的なFFTの処理実装の多くは、再帰的な処理を、系統だった再帰をしないアルゴリズムにより実現している。

クーリー–テューキー型アルゴリズムは変換をより小さい変換に分解していくので、後述のような他の離散フーリエ係数のアルゴリズムと任意に組み合わせることができる。とりわけ、//N// ≤ 8}} あたりまで分解すると、固定次数の高速なアルゴリズムに切り替えることが多い。

2.1.1. 原理の簡単な説明
---------------

離散フーリエ係数は、1}}の原始 N}} 乗根の1つ 1=//Wsub を使うと、次のように表せる。//例えば、1=//N// = 4}} のとき、、とすれば、離散フーリエ係数は行列を用いて表現すると

となる。入力列 xsub を添字の偶奇で分けて、以下のように変形する。



すると、サイズ 2}} のFFTの演算結果を用いて表現でき、サイズの分割ができる。

また、この分割手順を図にすると蝶のような図になることから、バタフライ演算とも呼ばれる。

バタフライ演算は、計算機上ではビット反転で実現される。DSPの中には、このバタフライ演算のプログラムを容易にするため、ビット反転アドレッシングを備えているものがある。

2.1.2. 原理の説明
------------

1=//N// = //PQ//}} とする。N}} 次離散フーリエ変換を、以下のようにP}} 次離散フーリエ変換とQ}} 次離散フーリエ変換に分解する。

N}} 次離散フーリエ変換：

を、1=//n// = 0, 1, ..., //N// − 1}} について計算することを考える。//n//}}, //k//}} を次のように書き換える。ただし 0 ≤ //n// ≤ //N// − 1}} また 0 ≤ //k// ≤ //N// − 1}} である。

すると

ここで、

と置くと、

となる。即ち、1=//F// = //F//}} の計算は、次の2ステップになる。

ステップ1: 1=//p// = 0, 1, ..., //P// − 1}} と 1=//r// = 0, 1, ..., //Q// − 1}} についてを計算する。これは、Q}}次の離散フーリエ変換の実行と、回転因子 exp}} の掛け算を、全ての p}}, r}} の組に対して行うことと見ることができる。ステップ2: 1=//s// = 0, 1, ..., //P// − 1}} と 1=//r// = 0, 1, ..., //Q// − 1}} についてを計算する。ここで、右辺は r}} を固定すれば、P}} 次の離散フーリエ変換である。ステップ1、2は、1=//N// = //PQ//}} 次の離散フーリエ変換を、Q}} 次の離散フーリエ変換と回転因子の掛け算の実行により、Q}} 組  の P}} 次離散フーリエ変換に分解したと見ることができる。

1=//N// = //Qk//}}  の場合には、上を繰り返せば、Q}} 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換にまで下げると、//F//}}を求めることができる。特に、Q}} が2}}または4}}の場合は、Q}}次の離散フーリエ変換は非常に簡単な計算になる。

1=//Q// = 2}} の場合は、exp}} は 1}} か −1}} なので、Q}} 次の離散フーリエ変換は符号の逆転と足し算だけで計算できる。
1=//Q// = 4}} の場合は、exp}} は 1}}, −1}}, i}}, −//i//}} のいずれかなので、Q}} 次の離散フーリエ変換の計算は、符号の逆転、実部虚部の交換と足し算だけで計算できる。このため、2}}の累乗あるいは4}}の累乗次の離散フーリエ変換は簡単に計算できる。実務的に用いられるのは、1=//Q// = 2}} か 1=//Q// = 4}} の場合のみである。なお、1=//Q// = 2}} か1=//Q// = 4}} の場合のこの部分のQ}}次の離散フーリエ変換のことを、バタフライ演算と言う。

また、1=//Q// = 2}}か1=//Q// = 4}}の場合において、計算を終了するまでに何回の「掛け算」が必要かを考える。符号の逆転、実部虚部の交換は「掛け算」として数えなければ、回転因子の掛け算のみが「掛け算」である。1=//N// = //Qsupの次数を1落とすためにN}}回の「掛け算」が必要であり、次数を//k//}}から0}}に落とすにはそれを//k//}}回繰り返す必要があるため、「掛け算」の数は 1=//Nk// = //N// logsub となる。高速フーリエ変換の計算において時間がかかるのは「掛け算」の部分であるため、これが「高速フーリエ変換では計算速度は //O//}} になる」ことの根拠になっている。

2.1.3. ビットの反転
-------------

上記の説明で、 の場合、1=//N// = //Qsup 個のデータから、1=//N// = //Qsup 個の計算結果

を計算する場合に、メモリの節約のため、0 ≤ //q// ≤ //Q// − 1}} と 0 ≤ //r// ≤ //Q// − 1}} を利用し、計算結果 を元データ のあった場所に格納することが多い。これが次の次数 //Q//sup でも繰り返されるため、とすると、次の次数の計算結果はのあった場所に格納される。繰り返せば、とすると、計算結果はのあった場所に格納される。

一方、

を、r}} を固定し s}} を変数とした //Q//sup 次離散フーリエ変換と見なして、とすると、

となる。繰り替えせば、

となるが、左辺について

より 1=//ssub, また右辺について//

より 1=//psub。このため、//

これは のあった場所に格納されている。

このように、求める解 が のあった場所に格納されていることを、ビット反転と言う。これは、Q}} 進法で表示した場合、 は となるのに対し、は逆から読んだとなるためである。

2.1.4. プログラムの例
--------------

以下は、高速フーリエ変換のプログラムを 1=//Q// = 4}} の場合にMicrosoft Visual Basicの文法を用いて書いた例である。

Const pi As Double = 3.14159265358979 '円周率 Dim Ndeg As Long '4^deg Dim Pdeg As Long '4^ Dim CR As Double '入力実数部 Dim CI As Double '入力虚数部 Dim FR As Double '出力実数部 Dim FI As Double '出力虚数部

deg=5 '任意に設定。5ならN=4^5=1024で計算 Ndeg=4^deg ReDim CR As Double '入力実数部 ReDim CI As Double '入力虚数部 ReDim FR As Double '出力実数部 ReDim FI As Double '出力虚数部 'ここで、変換される関数の実部をCRからCRに、虚部をCIからCIに入力しておくこと

'フーリエ変換 For i = 1 To deg Pdeg = 4 ^  For j0 = 0 To 4 ^  - 1 For j1 = 0 To Pdeg - 1 j = j1 + j0 * Pdeg * 4 'バタフライ演算 w1 = CR + CR + CR + CR w2 = CI + CI + CI + CI w3 = CR + CI - CR - CI w4 = CI - CR - CI + CR w5 = CR - CR + CR - CR w6 = CI - CI + CI - CI w7 = CR - CI - CR + CI w8 = CI + CR - CI - CR CR = w1 CI = w2 CR = w3 CI = w4 CR = w5 CI = w6 CR = w7 CI = w8 '回転因子 For k = 0 To 3 w1 = Cos w2 = -Sin w3 = CR * w1 - CI * w2 w4 = CR * w2 + CI * w1 CR(j + k
* Pdeg) = w3 CI = w4 Next k Next j1 Next j0 Next i 'ビット反転 For i = 0 To Ndeg - 1 k = i k1 = 0 For j = 1 To deg k1 = k1 +  * 4) * 4 ^  k = Int Next j FR = CR FI=CI Next i<WtXmlEndTag />

この例では、最深部 の繰り返し回数が Ndeg logsub Ndeg となっている。

2.2. その他のアルゴリズム
---------------

Prime Factor Algorithm 
レーダーのFFTアルゴリズム
Bluestein's FFT algorithm 　任意長のデータ列に対する変換が高速に可能である。
オドリツコ・ショーンハーゲ法 - アンドリュー・オドリツコ、アーノルド・ショーンハーゲ。
FFTW

3. 実数および対称的な入力への最適化
-------------------

多くの応用において、FFTに対する入力データは実数の列であり、このとき変換された出力の列は次の対称性を満たす：

そこで、多くの効率的なFFTアルゴリズム は入力データが実数であることを前提に設計されている。

入力データが実数の場合の効率化の手段としては、次のようなものがある。

クーリー-テューキー型アルゴリズムなど典型的なアルゴリズムを利用して、時間とメモリーの両方のコストを低減する。
入力データが偶数の長さのフーリエ係数はその半分の長さの複素フーリエ係数として表現できることを利用する。かつては実数の入力データに対するフーリエ係数を求めるのには、実数計算だけで行える離散ハートリー変換 を用いると効率的であろうと思われていた。しかしその後に、最適化された離散フーリエ変換  アルゴリズムの方が、離散ハートリー変換アルゴリズムに比べて必要な演算回数が少ないということが判明した。また当初は、実数入力に対してブルーン  FFT アルゴリズムは有利であると云われていたが、その後そうではないことが判った。

また、偶奇の対称性を持つ実入力の場合には、DFTはDCTや離散サイン変換となるので、演算と記憶に関してほぼ2倍の効率化が得られる。よって、そのような場合にはDFTのアルゴリズムをそのまま適用するよりも、DCTやDSTを適用してフーリエ係数を求める方が効率的である。

4. 応用
-----

スペクトラムアナライザ
OFDM変復調器OFDMの実装は、LSI化されたFFTおよびIFFTをそれぞれ復調器および変調器を用いて行われている。
フーリエ変換NMR核磁気共鳴  スペクトルを得るために使用される。
コンピュータ断層撮影 、核磁気共鳴画像法  等受像素子を360度回転させながら連続撮影した映像をフーリエ変換する事により、回転面の透過画像を合成する。
多倍精度の乗除算
自動列車停止装置 
FFTアナライザ周波数の分布を調べるために使用される。以前はハードウェアで信号を処理していたが、近年はCPUの性能が向上した為ソフトウェアで処理される。ノートパソコンとUSBで接続して使用するもの や、近年はデジタルオシロスコープにFFTの機能を内蔵している物もある。
電波天文学FX型デジタル分光相関器等を使用して星間分子のスペクトルを解析する。

5. 歴史
-----

高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー  とジョン・テューキー  が発見した とされているクーリー–テューキー型FFTアルゴリズムを呼ぶ。同時期に高橋秀俊がクーリーとテューキーとは全く独立にフーリエ変換を高速で行うためのアルゴリズムを考案していた。しかし、1805年頃に既にガウスが同様のアルゴリズムを独自に発見していた。ガウスの論文以降、地球物理学や気候や潮位解析などの分野などで測定値に対する調和解析は行われていたので、計算上の工夫を必要とする応用分野で受け継がれていたようである。 以下の書籍にも、天体観測の軌道の補間のためにガウスが高速フーリエ変換を利用したことが書かれている。

Elena Prestini:"The Evolution of Applied Harmonic Analysis", Springer, ISBN 978-0-8176-4125-2 のSec.3.10 'Gauss and the asteroids: history of the FFT'.

6. 参考文献
-------

7. 関連記事
-------

フーリエ変換
離散フーリエ変換 
Overlap–add method /
スペクトラムアナライザ
FFTPACK
時系列
Math Kernel Library
ショーンハーゲ・ストラッセン法
Sparse Fourier transform ※ スパース・フーリエ変換

8. 学習用図書
--------

今後記述を追加の予定

Henri J. Nussbaumer: "Fast Fourier Transform and Convolution Algorithms",2nd Ed.,Springer-Verlag, ISBN 978-3-540-11825-1 .
E.Oran Brigham：「高速フーリエ変換」、科学技術出版社 .
Henri J. Nussbaumer:「高速フーリエ変換のアルゴリズム」、科学技術出版社、ISBN 978-4876530069 ．
William L. Briggs and Van Emden Henson: "The DFT: An Owners' Manual for the Discrete Fourier Transform", SIAM, ISBN 978-0-898713-42-8 .
Gerlind Plonka, Daniel Potts, Gabriele Steidl and Manfred Tasche: "Numerical Fourier Analysis", Birkhaeuser, ISBN 978-3030043056 (2019年2月）