高速フーリエ変換

提供: miniwiki
移動先:案内、検索

高速フーリエ変換（こうそくフーリエへんかん、英: fast Fourier transform, FFT）は、離散フーリエ変換（英: discrete Fourier transform, DFT）を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換（英: inverse fast Fourier transform, IFFT）と呼ぶ。

歴史

高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー（English版） (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した[1] とされているクーリー–テューキー型FFTアルゴリズム（English版）を呼ぶ[2]。同時期に高橋秀俊がクーリーとテューキーとは全く独立にフーリエ変換を高速で行うためのアルゴリズムを考案していた[3]。しかし、1805年頃に既にガウスが同様のアルゴリズムを独自に発見していた[4]。

概要

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

[math]F(t)= \sum_{x=0}^{N-1} f(x) e^{-i\frac{2 \pi t x}{N}} .[/math]

このとき、{x = 0, 1, 2, ..., N − 1} を標本点と言う。

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

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

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

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

逆変換

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

離散フーリエ変換を

[math]F(t) = \sum_{x=0}^{N-1} f(x) e^{-i\frac{2 \pi t x}{N}}[/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} }}[/math]

となる。

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

(1) 複素共役を取り、F(t) を求める、
(2) F(t) の正変換の離散フーリエ変換を高速フーリエ変換で行う、
(3) その結果の複素共役を取り、N で割る

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

アルゴリズム

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

"「クーリー–テューキー型FFTアルゴリズム（English版）」"

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

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

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

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

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

原理の簡単な説明

ファイル:FFT visual 12.jpg
データ数1２の離散フーリエ変換の模式図。時計を模した図形は1の12乗根の一つを表している。時計の針の向きと色は1の12乗根の偏角を表す。この図で表される行列をデータ列にかけることで離散フーリエ変換が得られる。上図で表されるような列の並び替えを行うことで、元の行列のパターンはデータ数6の離散フーリエ変換のパターンに分解できる。この繰り返しにより最終的にはデータ数３のフーリエ変換に帰着される。
ファイル:FFT visual 100.jpg
データ数100の離散フーリエ変換の模式図。色は1の100乗根の偏角を表す。バタフライ演算により元の行列のパターンは最終的にデータ数5の離散フーリエ変換のパターンに分解される。
ファイル:FFT Butterfly radix8.svg
FFTのバタフライ演算

離散フーリエ係数は、1の原始 N 乗根の1つ WN = e−2πi/N を使うと、次のように表せる。

[math]F(t) = \sum_{x=0}^{N-1} f(x) W_N^{tx} .[/math]

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

[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]

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

[math] \begin{align} \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^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_1\\X_2\\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]

となる。即ち、F(n) = F(sQ + r) の計算は、次の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]
の実行と、回転因子 exp(−2πipr/PQ) の掛け算を、全ての 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 = Qk (P = Qk − 1) の場合には、上を繰り返せば、Q 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換（何もしないことと同じ）にまで下げると、F(t)を求めることができる。特に、Q が2または4の場合は、Q次の離散フーリエ変換は非常に簡単な計算になる。

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

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

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

ビットの反転

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

[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 を利用し、計算結果 [math]f_1(p,r)[/math] を元データ[math]f(rQ^{k-1}+p)[/math] のあった場所に格納することが多い。これが次の次数 Qk − 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 を変数とした Qk − 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\lt Q^k[/math]

より sk = 0, また右辺について

[math]Q^{k-k}-1=0[/math]

より pk = 0。このため、

[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]となるためである。

プログラムの例

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

Const pi As Double = 3.14159265358979   '円周率
Dim Ndeg As Long '4^deg
Dim Pdeg As Long '4^(deg-i)
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(Ndeg - 1) As Double '入力実数部
ReDim CI(Ndeg - 1) As Double '入力虚数部
ReDim FR(Ndeg - 1) As Double '出力実数部
ReDim FI(Ndeg - 1) As Double '出力虚数部
'ここで、変換される関数の実部をCR(0)からCR(Ndeg-1)に、虚部をCI(0)からCI(Ndeg-1)に入力しておくこと

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

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

その他のアルゴリズム


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

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

[math]F(-t) = \overline{F(t)} .[/math]

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

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

  • クーリー–テューキー型アルゴリズムなど典型的なアルゴリズムを利用して、時間とメモリーの両方のコストを低減する。
  • 入力データが偶数の長さのフーリエ係数はその半分の長さの複素フーリエ係数として表現できる（出力の実数/虚数成分は、それぞれ入力の偶関数/奇関数成分に対応する）ことを利用する。

かつては実数の入力データに対するフーリエ係数の計算は離散ハートリー変換（English版） (discrete Hartley transform, DHT) を用いればさらに効率化できると思われていた。しかしその後に、最適化された離散フーリエ変換 (discrete Fourier transform, DFT) アルゴリズムの方が、離散ハートリー変換アルゴリズムに比べて必要な演算回数が少ないことが判明した。また、ブルーン (Bruun) FFT アルゴリズムは実数入力に対して有利であると最初は云われていたが現在ではそうではない。

また、偶奇の対称性を持つ実入力の場合にはさらに最適化ができて、DFTはDCTやDST（English版）となり、時間とメモリーに関してほぼ2倍の高速化が得られるから、そのような場合にはDFTのアルゴリズムを適用するのではなく、DCTやDSTを直接適用した計算によりフーリエ係数を求める方が良い。

応用

OFDM（日本および欧州で地上デジタルテレビジョン放送やADSL等に用いられる変調方式）の実装は、LSI化されたFFTおよびIFFT（逆変換）をそれぞれ復調器および変調器を用いて行われている。
核磁気共鳴 (NMR) スペクトルを得るために使用される。
受像素子を360度回転させながら連続撮影した映像をフーリエ変換する事により、回転面の透過画像を合成する。
周波数の分布を調べるために使用される。以前はハードウェアで信号を処理していたが、近年はCPUの性能が向上した為ソフトウェアで処理される。ノートパソコンとUSBで接続して使用するもの[6] や、近年はデジタルオシロスコープにFFTの機能を内蔵している物もある。
FX型デジタル分光相関器等を使用して星間分子のスペクトルを解析する[7][8]。

参考文献

  1. ↑ 1.0 1.1 J. W. Cooley and J. W. Tukey: Math. of Comput. 19 (1965) 297.
  2. ↑ IEEE Archives: History of FFT with Cooley and Tukey.
  3. ↑ 『東京大学大型計算機センターニュース』第2巻Supplement 2、1970年。
  4. ↑ Carl Friedrich Gauss, "Nachlass: Theoria interpolationis methodo nova tractata", Werke band 3, 265–327 (Konigliche Gesellschaft der Wissenschaften, Gottingen, 1866). See also M. T. Heideman, D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform", IEEE ASSP Magazine 1 (4), 14–21 (1984).
  5. ↑ 例えば、H. V. Sorensen, D. L. Jones, M. T. Heideman, and C. S. Burrus, "Real-valued fast Fourier transform algorithms," IEEE Trans. Acoust. Speech Sig. Processing ASSP-35, 849–863 (1987).
  6. ↑ FFT spectrum analyzer
  7. ↑ 惑星大気の観測「SPART」
  8. ↑ 空間FFT電波干渉計による電波天体の高速撮像

関連記事

外部リンク



「http:///http:///mymemo.xyz/mymemo.xyz/wiki/index.php?title=高速フーリエ変換&oldid=114640」から取得