高速フーリエ変換

高速フーリエ変換（こうそくフーリエへんかん、テンプレート:Lang-en-short、FFT）とは、離散フーリエ変換 (Discrete Fourier Transform、DFT) を計算機上で高速に計算するアルゴリズム。逆変換をIFFT (Inverse FFT) という。

歴史

高速フーリエ変換といえば一般的には1965年、テンプレート:仮リンク (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した^[1]とされているテンプレート:仮リンクを呼ぶ^[2]。しかし、1805年ごろにガウスが同様のアルゴリズムを独自に発見していた^[3]。

概要

複素関数f(x)の離散フーリエ変換 (DFT)である複素関数F(t) は以下で定義される。テンプレート:Indent \quad \quad </math>}} このとき、<math>{x=0,1,2,\dots, N-1}</math>を標本点と言う。

これを直接計算したときの時間計算量は O (N² ) である（Oはランダウの記号）。

高速フーリエ変換は、この結果を、次数 N が2の累乗のときに O (N log N ) の計算量で得るアルゴリズムである。より一般的には、次数が <math>N = \prod n_i</math> と素因数分解できるとき、<math>O ( N \sum n_i ) </math> の計算量となる。次数が2の累乗のときが最も高速に計算でき、アルゴリズムも単純になるので、0詰めで次数を調整することもある。

高速フーリエ変換を使って、畳み込み積分などの計算を高速に求めることができる。これも計算量をO (N² ) → O (N log N ) に落とせる。

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

逆変換

逆変換 (IFFT) は正変換 (FFT) と同じと考えて良いが、指数の符号が逆であり、係数<math>\frac{1}{N}</math>が掛かる。

離散フーリエ変換を以下で定義したとき、テンプレート:Indent \quad \quad </math>}} 逆変換は、

<math> f(x) = \frac{1}{N} \sum_{t=0}^{N-1} F(t) e^{i\frac{2\pi{tx}}{N}}= \frac{1}{N} \overline{ \sum_{t=0}^{N-1} \overline{F(t)} e^{-i\frac{2 \pi t x}{N} }} \quad \quad</math>

となる。

このため、<math> F(t) </math>の離散フーリエ逆変換を求めるには、

(1)複素共役を取り、<math> \overline{F(t)} </math>を求める、

(2)<math> \overline{F(t)} </math>の正変換の離散フーリエ変換を高速フーリエ変換で行う、

(3)その結果の複素共役を取り、Nで割る

とすれば良く、正変換の高速フーリエ変換のプログラムがあれば、逆変換は容易に作ることができる。

アルゴリズム

Cooley-Tukey型FFTアルゴリズム

テンプレート:Main Cooley-Tukey型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。

分割統治法 (divide and conquer) を使ったアルゴリズムで、 N = N₁ N₂ のサイズの変換を、より小さいサイズである N₁ 、N₂ のサイズの変換に分割していくことで高速化を図っている。

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

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

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

原理の簡単な説明

ファイル:FFT Butterfly radix8.svg

FFTのバタフライ演算

離散フーリエ係数は、1の原始 N 乗根の1つ <math>W_{N} \equiv e^{-i\frac{2\pi}{N}}</math>を使うと、次のように表せる。

<math>F(t) = \sum_{x=0}^{N-1} f(x) W_N^{tx} \qquad </math>

例えば、N = 4 のとき、離散フーリエ係数は行列を用いて表現すると（W ≡ W₄ と略記）

<math>

\begin{bmatrix}X_0\\X_1\\X_2\\X_3\end{bmatrix}= \begin{bmatrix} W^0&W^0&W^0&W^0\\ W^0&W^1&W^2&W^3\\ W^0&W^2&W^4&W^6\\ W^0&W^3&W^6&W^9 \end{bmatrix} \begin{bmatrix}x_0\\x_1\\x_2\\x_3\end{bmatrix} </math> となる。入力列 x_k を添字の偶奇で分けて、以下のように変形する。

<math>

\begin{align} \begin{bmatrix}X_0\\X_2\\X_1\\X_3\end{bmatrix} & = \begin{bmatrix} W^0&W^0&W^0&W^0\\ W^0&W^2&W^1&W^3\\ W^0&W^4&W^2&W^6\\ W^0&W^6&W^3&W^9 \end{bmatrix} \begin{bmatrix}x_0\\x_2\\x_1\\x_3\end{bmatrix} \\ & = \begin{bmatrix} W^0&W^0&W^0W^0&W^0W^0\\ W^0&W^2&W^1W^0&W^1W^2\\ W^0&W^0&W^2W^0&W^2W^0\\ W^0&W^2&W^3W^0&W^3W^2 \end{bmatrix} \begin{bmatrix}x_0\\x_2\\x_1\\x_3\end{bmatrix} \end{align} </math> すると、サイズ2のFFTの演算結果を用いて表現でき、サイズの分割ができる。

<math>\begin{bmatrix}X_0\\X_2\\X_1\\X_3\end{bmatrix}=

\begin{bmatrix} 1&0&W^0&0\\ 0&1&0&W^1\\ 1&0&W^2&0\\ 0&1&0&W^3 \end{bmatrix}\, \begin{bmatrix} W_2^0&W_2^0&0&0\\ W_2^0&W_2^1&0&0\\ 0&0&W_2^0&W_2^0\\ 0&0&W_2^0&W_2^1 \end{bmatrix}\, \begin{bmatrix} x_0\\ x_2\\ x_1\\ x_3 \end{bmatrix} </math> また、この分割手順を図にすると蝶のような図になることから、バタフライ演算とも呼ばれる。

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

原理の説明

N =PQ とする。N 次離散フーリエ変換：

<math>F(n) = \sum_{k=0}^{N-1}f(k)\exp\left(-2\pi i\frac{nk}{N}\right)</math>

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

<math>\begin{align}

n &= sQ+r \qquad 0\le s\le P-1,\,0\le r\le Q-1 \\ k &= qP+p \qquad 0\le p\le P-1,\;0\le q\le Q-1 \end{align}</math> すると

<math>\begin{align}

F(n) &= F(sQ+r) = \sum_{p=0}^{P-1}\sum_{q=0}^{Q-1}f(qP+p)\exp\left(-2\pi i\frac{(qP+p)(sQ+r)}{PQ}\right) \\ &= \sum_{p=0}^{P-1}\sum_{q=0}^{Q-1}f(qP+p)\exp\left(-2\pi i\left(\frac{rq}{Q}+\frac{sp}{P}+\frac{pr}{PQ}\right)\right) \\ &= \sum_{p=0}^{P-1}\left[\exp\left(-2\pi i\left(\frac{sp}{P}+\frac{pr}{PQ}\right)\right)\sum_{q=0}^{Q-1}f(qP+p)\exp\left(-2\pi i\frac{rq}{Q}\right)\right]\end{align}</math>

ここで、

<math>f_1(p,r) = \exp\left(-2\pi i\frac{pr}{PQ}\right)\sum_{q=0}^{Q-1}f(qP+p)\exp\left(-2\pi i\frac{rq}{Q}\right)</math>

と置くと、

<math>F(n) = F(sQ+r) = \sum_{p=0}^{P-1}\exp\left(-2\pi i\frac{sp}{P}\right)f_1(p,r)</math>

となる。即ち、<math>F(n)=F(sQ+r)</math>の計算は、次の2ステップになる。

ステップ1

p = 0, 1, ... , P -1、r = 0, 1, ... , Q -1 について

<math>f_1(p,r)= \exp\left(-2\pi i\frac{pr}{PQ}\right)\sum_{q=0}^{Q-1}f(qP+p)\exp\left(-2\pi i\frac{rq}{Q}\right)</math>

を計算する。これは、Q 次の離散フーリエ変換

<math>\sum_{q=0}^{Q-1}f(qP+p)\exp\left(-2\pi i\frac{rq}{Q}\right)</math>

の実行と、回転因子<math>\exp(-2\pi ipr/PQ)</math>の掛け算を、全てのp 、r の組（PQ = N 通り）に対して行うことと見ることができる。

ステップ2

s = 0, 1, ... , P -1、r = 0, 1, ... , Q -1 について

<math>F(sQ+r) = \sum_{p=0}^{P-1}\exp\left(-2\pi i\frac{sp}{P}\right)f_1(p,r)</math>

を計算する。ここで、右辺はr を固定すれば、P 次の離散フーリエ変換である。

ステップ1、2は、N = PQ 次の離散フーリエ変換を、Q 次の離散フーリエ変換と回転因子の掛け算の実行により、Q 組（r = 0, 1, ... , Q -1）のP 次離散フーリエ変換に分解したと見ることができる。

N = Q^k （P = Q^{k -1} ）の場合には、上を繰り返せば、Q 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換（何もしないことと同じ）にまで下げると、F (t )を求めることができる。特に、Q が2または4の場合は、Q 次の離散フーリエ変換は非常に簡単な計算になる。

Q = 2の場合は、<math>\exp(-2\pi i\,rq/Q)</math>は1か-1なので、Q 次の離散フーリエ変換は符号の逆転と足し算だけで計算できる。
Q = 4の場合は、<math>\exp(-2\pi i\,rq/Q)</math>は1, -1, i, -iのいずれかなので、Q 次の離散フーリエ変換の計算は、符号の逆転、実部虚部の交換と足し算だけで計算できる。

このため、2の累乗あるいは4の累乗次の離散フーリエ変換は簡単に計算できる。実務的に用いられるのは、Q = 2かQ = 4の場合のみである。なお、Q = 2かQ = 4の場合のこの部分のQ 次の離散フーリエ変換のことを、バタフライ演算と言う。

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

ビットの反転

上記の説明で、N = Q^k （P = Q^{k -1} ）の場合、N = Q^k 個のデータ<math>f(qQ^{k-1}+p)</math>から、N = Q^k 個の計算結果

<math>f_1(p,r)= \exp\left(-2\pi i\frac{pr}{Q^k}\right)\sum_{q=0}^{Q-1}f(qQ^{k-1}+p)\exp\left(-2\pi i\frac{rq}{Q}\right)</math>

を計算する場合に、メモリの節約のため、0≤q ≤Q -1、0≤r ≤Q -1 を利用し、計算結果f₁(p ,r ) を元データ<math>f(rQ^{k-1}+p)</math> のあった場所に格納することが多い。これが次の次数Q^{k -1} でも繰り返されるため、<math>p=q_2Q^{k-2}+p_2</math>とすると、次の次数の計算結果<math>f_2(p_2,q_2,q)</math>は<math>f(qQ^{k-1}+q_2Q^{k-2}+p_2)</math>のあった場所に格納される。繰り返せば、<math>t=q_1Q^{k-1}+q_2Q^{k-2}+\cdots+q_k</math>とすると、計算結果<math>f_k(p_k,q_k,q_{k-1},\dots,q_2,q_1)</math>は<math>f(q_1Q^{k-1}+q_2Q^{k-2}+\cdots+q_{k-1}Q+p_k)</math>のあった場所に格納される。

一方、

<math>F(sQ+r)=\sum_{p=0}^{Q^{k-1}-1}\exp\left(-2\pi i\frac{sp}{Q^{k-1}}\right)f_1(p,r)</math>

を、r を固定しs を変数としたQ^{k -1} 次離散フーリエ変換と見なして、<math>s=s_2Q+r_2</math>とすると、

<math>F(s_2Q^2+r_2Q+r) = \sum_{p_2=0}^{Q^{k-2}-1}\exp\left(-2\pi i\frac{s_2p_2}{Q^{k-2}}\right)f_2(p_2,r_2,r)</math>

となる。繰り替えせば、

<math>F(s_kQ^k+r_kQ^{k-1}+\cdots+r_2Q+r_1) = \sum_{p_k=0}^{Q^{k-k}-1}\exp\left(-2\pi i\frac{s_kp_k}{Q^{k-k}}\right)f_k(p_k,r_k,r_{k-1},\dots,r_2,r_1)</math>

となるが、左辺について

<math>s_kQ^k+r_kQ^{k-1}+\cdots+r_2Q+r_1<Q^k</math>

より<math>s_k=0</math>、また右辺について

より<math>p_k=0</math>。このため、

<math>F(r_kQ^{k-1}+\cdots+r_2Q+r_1)=f_k(0,r_k,r_{k-1},\dots,r_2,r_1)</math>

これは<math>f(r_1Q^{k-1}+r_2Q^{k-2}+\cdots+r_{k-1}Q+r_k)</math>のあった場所に格納されている。

このように、求める解<math>F(r_kQ^{k-1}+\cdots+r_2Q+r_1)</math> が<math>f(r_1Q^{k-1}+r_2Q^{k-2}+\cdots+r_{k-1}Q+r_k)</math> のあった場所に格納されていることを、ビット反転と言う。これは、Q 進数で表示した場合、<math>r_kQ^{k-1}+\cdots+r_2Q+r_1</math> は <math>(r_kr_{k-1}\dots r_2r_1)_Q</math>となるのに対し、<math>r_1Q^{k-1}+r_2Q^{k-2}+\cdots+r_{k-1}+r_k</math>は逆から読んだ<math>(r_1r_2\dots r_{k-1}r_k)_Q</math>となるためである。

プログラムの例

Microsoft Visual Basicの文法で、高速フーリエ変換のプログラムを書くと、次のような例が考えられる。この例ではQ = 4である。

（プログラムはじまり）</br>
Const pi As Double = 3.14159265358979   '円周率</br>
Dim Ndeg As Long '4^deg</br>
Dim Pdeg As Long '4^(deg-i)</br>
Dim CR() As Double   '入力実数部</br>
Dim CI() As Double   '入力虚数部</br>
Dim FR() As Double   '出力実数部</br>
Dim FI() As Double   '出力虚数部</br>
</br>
deg=5 '任意に設定。5ならN=4^5=1024で計算</br>
Ndeg=4^deg</br>
ReDim CR(Ndeg - 1) As Double '入力実数部</br>
ReDim CI(Ndeg - 1) As Double '入力虚数部</br>
ReDim FR(Ndeg - 1) As Double '出力実数部</br>
ReDim FI(Ndeg - 1) As Double '出力虚数部</br>
'ここで、変換される関数の実部をCR(0)からCR(Ndeg-1)に、虚部をCI(0)からCI(Ndeg-1)に入力しておくこと</br>
</br >
'フーリエ変換</br>
For i = 1 To deg</br>
 Pdeg = 4 ^ (deg - i)</br>
 For j0 = 0 To 4 ^ (i - 1) - 1</br>
  For j1 = 0 To Pdeg - 1</br>
   j = j1 + j0 * Pdeg * 4</br>
   'バタフライ演算(Q=4)</br>
   w1 = CR(j) + CR(j + Pdeg) + CR(j + 2 * Pdeg) + CR(j + 3 * Pdeg)</br>
   w2 = CI(j) + CI(j + Pdeg) + CI(j + 2 * Pdeg) + CI(j + 3 * Pdeg)</br>
   w3 = CR(j) + CI(j + Pdeg) - CR(j + 2 * Pdeg) - CI(j + 3 * Pdeg)</br>
   w4 = CI(j) - CR(j + Pdeg) - CI(j + 2 * Pdeg) + CR(j + 3 * Pdeg)</br>
   w5 = CR(j) - CR(j + Pdeg) + CR(j + 2 * Pdeg) - CR(j + 3 * Pdeg)</br>
   w6 = CI(j) - CI(j + Pdeg) + CI(j + 2 * Pdeg) - CI(j + 3 * Pdeg)</br>
   w7 = CR(j) - CI(j + Pdeg) - CR(j + 2 * Pdeg) + CI(j + 3 * Pdeg)</br>
   w8 = CI(j) + CR(j + Pdeg) - CI(j + 2 * Pdeg) - CR(j + 3 * Pdeg)</br>
   CR(j) = w1</br>
   CI(j) = w2</br>
   CR(j + Pdeg) = w3</br>
   CI(j + Pdeg) = w4</br>
   CR(j + 2 * Pdeg) = w5</br>
   CI(j + 2 * Pdeg) = w6</br>
   CR(j + 3 * Pdeg) = w7</br>
   CI(j + 3 * Pdeg) = w8</br>
   '回転因子</br>
   For k = 0 To 3</br>
    w1 = Cos(2 * pi * j * k / Pdeg / 4)</br>
    w2 = -Sin(2 * pi * j * k / Pdeg / 4)</br>
    w3 = CR(j + k * Pdeg) * w1 - CI(j + k * Pdeg) * w2</br>
    w4 = CR(j + k * Pdeg) * w2 + CI(j + k * Pdeg) * w1</br>
    CR(j + k * Pdeg) = w3</br>
    CI(j + k * Pdeg) = w4</br>
   Next k</br>
  Next j1</br>
 Next j0</br>
Next i</br>
'ビット反転</br>
For i = 0 To Ndeg - 1</br>
 k = i</br>
 k1 = 0</br>
 For j = 1 To deg</br>
  k1 = k1 + (k - Int(k / 4) * 4) * 4 ^ (deg - j)</br>
  k = Int(k / 4)</br>
 Next j</br>
 FR(i) = CR(k1)</br>
 FI(i)=CI(k1)</br>
Next i</br>
（プログラムおわり）

この例では、最深部(For k、Next kの間の部分)の繰り返し回数がNdeg<math>\log_4</math>Ndegとなっている。

その他のアルゴリズム

テンプレート:節stub

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

多くの実用面において、FFTへの入力は実数のみの列(実入力)であり、このとき出力の列は次の対称性を満たす（*は複素共役）：

<math>F(-t) = \overline{F(t)} \quad </math>．

そして、多くの効率的なFFTアルゴリズム（例えば^[4]）はこの条件を前堤に設計されている。

実入力への効率化には、次のような手段がある。

Cooley-Tukey型アルゴリズムなど典型的なアルゴリズムを利用して、時間とメモリーの両方の浪費を低減する。
偶数区分数のフーリエ係数はその半分の長さの複素フーリエ係数として表現することができる（出力の実数/虚数成分は、それぞれ入力の偶関数/奇関数成分に対応する）ことを利用する。

かつては実入力に対するフーリエ係数はテンプレート:仮リンク (DHT : Discrete Hartley transform) でより効率的に計算できると思われていた。しかし、その後最適化された離散フーリエ変換(DFT：Discrete Fourier Transform)アルゴリズムが、離散ハートリー変換アルゴリズムよりも必要とする演算回数を下まわることがわかった。また、Bruun's FFTアルゴリズムは、はじめこそ実入力に対して有利といわれていたが、今ではそういわれてはいない。

さらに偶奇の対称性を持つ実入力の場合にはさらに最適化ができ、DFTはDCTやテンプレート:仮リンクになり、時間とメモリーの（ほとんど）2つに関して最適化が得られる。この場合には直接DFTのアルゴリズムを応用しないでも、DCTやDSTの計算によってフーリエ係数を求めることができる。

応用

スペクトラムアナライザ
OFDM変復調器

日本および欧州で地上デジタルテレビジョン放送に用いられる変調方式であるOFDMの実装は、LSI化されたFFTおよびIFFT（逆変換）をそれぞれ復調器および変調器を用いて行われている。

CT、MRI等

受像素子を360度回転させながら連続撮影した映像をフーリエ変換する事により、回転面の透過画像を合成する。

多倍精度の乗除算
自動列車停止装置 (例: JR西日本の最新型車両。地上子が発振する周波数に高速フーリエ変換を用いている)
FFTアナライザ

周波数の分布を調べるために使用される。近年はデジタルオシロスコープにFFTの機能を内蔵している物もある。また、以前はハードウェアで信号を処理していたが近年はCPUの性能が向上した為ソフトウェアで処理している。また、ノートPCとUSBで接続して使用するものもある。

参考文献

テンプレート:脚注ヘルプテンプレート:Reflist