Pythonプログラミング（ステップ８・関数・高速フーリエ変換）

このページでは、信号や画像の処理に欠かせない、高速フーリエ変換（FFT）のアルゴリズムについて考える。

１．離散フーリエ変換の再帰処理

$N$点のデータ $x_0, x_1, \cdots, x_{N-1}$ が与えられているとき、その離散フーリエ変換を $$ X_k = \sum_{n=0}^{N-1} x_n \exp\left( - \frac{2 \pi n k}{N} i \right) \tag{1} $$ で定義する（ここで、$k=0, 1, \cdots, N-1$。$i$は虚数単位）。 つまり、$N$点の系列データ$\{x_n\}$ を、「見方」を変えて、同じく$N$点の$\{X_k\}$に変換するわけである。

指数関数の部分は、オイラーの公式を思い出すと、 $$ \begin{eqnarray} \exp\left( - \frac{2 \pi n k}{N} i \right) & = & \cos\left( - \frac{2 \pi n k}{N} \right) + i \sin\left( - \frac{2 \pi n k}{N} \right) \\ & = & \cos\left(\frac{2 \pi n k}{N} \right) - i \sin\left(\frac{2 \pi n k}{N} \right) \end{eqnarray} $$ であるから、この変換は、信号にコサインとサインの波を掛け合わせながら総和を取ることで、信号の中から角振動数 $2 \pi k / N$ の成分を抽出する操作であることを こちらのページですでに述べた。

ここで、データ点数$N$が2のべき乗であるような場合（$N=2,4,8,16, \cdots$)について考えてみる。

最も単純な場合として、１点からなるデータ $x_0$ のフーリエ変換は $X_0 = x_0$、２点からなるデータ $x_0, x_1$ のフーリエ変換は $$ \begin{eqnarray} X_0 & = & x_0 \exp\left( - \frac{2 \pi i \times 0 \times 0}{2} \right) + x_1 \exp\left( - \frac{2 \pi i \times 0 \times 1}{2} \right) = x_0 + x_1 \\ X_1 & = & x_0 \exp\left( - \frac{2 \pi i \times 1 \times 0}{2} \right) + x_1 \exp\left( - \frac{2 \pi i \times 1 \times 1}{2} \right) = x_0 - x_1 \end{eqnarray} $$ である。

データ点数がさらに多い場合、(1)式の総和を、偶数番目と奇数番目の項に分けると $$ \begin{eqnarray} X_k & = & \sum_{m=0}^{N/2-1} x_{2m} \exp\left( - \frac{2 \pi (2m) k}{N} i \right) + \sum_{m=0}^{N/2-1} x_{2m+1} \exp\left( - \frac{2 \pi (2m+1) k}{N} i \right) \\ & = & \sum_{m=0}^{N/2-1} x_{2m} \exp\left( - \frac{2 \pi m k}{N/2} i \right) + \sum_{m=0}^{N/2-1} x_{2m+1} \exp\left( - \frac{2 \pi m k}{N/2} i \right) \exp\left(- \frac{2\pi k}{N} i \right) \end{eqnarray} \tag{2} $$ となるが、 $$ X^{(e)}_k = \sum_{m=0}^{N/2-1} x_{2m} \exp\left( - \frac{2 \pi m k}{N/2} i \right) \tag{3} $$ は、$x_n$の偶数番目のみを抽出したデータ点数がちょうど半分のデータについての、同様に $$ X^{(o)}_k = \sum_{m=0}^{N/2-1} x_{2m+1} \exp\left( - \frac{2 \pi m k}{N/2} i \right) \tag{4} $$ のほうは$x_n$の奇数番目のみのデータのフーリエ変換にほかならない。 つまり、$k \lt N/2$ に対しては、$X_k$ を「半分ずつ」のデータのフーリエ変換の和 $$ X_k = X^{(e)}_k + X^{(o)}_k \exp\left(- \frac{2\pi k}{N} i \right) \tag{5} $$ に分解できる。

さらに、 $$ X_{N-k} = \sum_{n=0}^{N-1} x_n \exp\left( - \frac{2 \pi n (N-k)}{N} i \right) = \sum_{n=0}^{N-1} x_n \exp\left(\frac{2 \pi i n k}{N} i \right) = X^*_{k} \tag{6} $$ であるから、$k \lt N/2$について$X_k$を求めれれば、残りの部分は複素共役を取るだけで済む。

あるいは $$ \begin{eqnarray} X_{k+N/2} & = & \sum_{n=0}^{N-1} x_n \exp\left( - \frac{2 \pi n (k + N/2)}{N} i \right) \\ & = & \sum_{n=0}^{N-1} x_n \exp\left( - \frac{2 \pi n k}{N} i \right) \exp\left( - n \pi i \right) \\ & = & X^{(e)}_k - X^{(o)}_k \exp\left(- \frac{2\pi k}{N} i \right) \end{eqnarray} \tag{7} $$ から、$k \ge N/2$についての$X_k$を求めてもよい。

以上で述べた、計算過程に現れる再帰的なパターンと、前半と後半が複素共役となる関係を考慮して、 フーリエ変換を再帰的な関数でコーディングした例を以下に示す。

# coding: utf-8

import math
import cmath

def fft(x):
    N=len(x)
    if N==1:
        return [ x[0] ]
    else:
        XE = fft(x[0::2])
        XO = fft(x[1::2])
        X = [0]*N
        X[0] = XE[0] + XO[0]
        for k in range(1,N//2):
            X[k] = XE[k] + XO[k] * cmath.exp(-2*math.pi*k*(0+1j)/N)
            X[N-k] = X[k].conjugate()
        X[N//2] = XE[0] - XO[0] 
        return X

# メイン部 
N=64
x = [math.sin(2*math.pi*3*k/N)+math.cos(2*math.pi*5*k/N) for k in range(N) ]
X = fft(x)
print(X)

　練習：逆変換の実装

(1)式の離散フーリエ変換の逆変換は $$ x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \exp\left(\frac{2 \pi n k}{N} i \right) \tag{8} $$ で定義される。

例題コードをひな形にして、この逆変換を行う関数を設計してみなさい。

　練習：基数の変更

例題のコードをもとに、データ点数$N$が3のべき乗（$3,9,27,81,\cdots)$の場合に対応できるように拡張してみなさい。

２．ボトムアップな処理

次に、同じ計算を、少し異なる視点からもみておこう。

式(5)および(7)のルールに従って、$X_k$の計算をより小さなサイズのフーリエ変換に還元する過程を、$N=8$の場合について図示したのが以下である。

最終的にはサイズ1のフーリエ変換にたどり着くことになるが、その際の $x_n$ と $X^{(\cdots)}_0$ とには $$ \begin{eqnarray} x_0 \rightarrow X^{(eee)} \\ x_4 \rightarrow X^{(eeo)} \\ x_2 \rightarrow X^{(eoe)} \\ x_6 \rightarrow X^{(eoo)} \\ x_1 \rightarrow X^{(oee)} \\ x_5 \rightarrow X^{(oeo)} \\ x_3 \rightarrow X^{(ooe)} \\ x_7 \rightarrow X^{(ooo)} \\ \end{eqnarray} $$ という対応があることが分かる。 これは $x_n$ の添字 $n$ を２進数で表記し、0をeに、1をoに置き換えてから、その順序を反転させたパターンになっており、 2のべき乗であるような$N$については常に成り立つ。

そこで、上記の処理の流れに沿って、再帰を使わずボトムアップにフーリエ変換を行うコードの例の以下に示す：

# coding: utf-8

import math
import cmath

def bit_reverse(N,n):
    m=0
    while N>1:
        m *= 2
        m += n%2
        n = n//2
        N = N//2
    return m

def fft(x):
    N = len(x)
    W = cmath.exp(-2*math.pi*(0+1j)/N)
    X=[0]*N
    X2=[0]*N # work space
    # very initial loop
    for j in range(0,N):
        X[j]   = x[bit_reverse(N,j)]
    p=N//2
    M=2
    while M<=N:
        for k in range(0,N,M):
            for j in range(0,M//2):
                X2[k+j]      = X[k+j] + W**(p*j) * X[k+j+M//2]
                X2[k+j+M//2] = X[k+j] - W**(p*j) * X[k+j+M//2]
            for j in range(0,M):
                X[k+j] = X2[k+j]
        M = M*2
        p = p//2
    return X

# メイン部
N=64
x = [math.sin(2*math.pi*3*k/N)+math.cos(2*math.pi*5*k/N) for k in range(N) ]
X = fft(x)
print(X)

　練習：インプレースな処理

上記のコードでは、関数fft()で、入力データを格納する配列x、変換結果を格納する配列 Xに加え、 途中計算結果を保持するための配列X2を用いている。

さらに使用するメモリの量を削減する余地はないか、検討してしてみなさい。