FFTの原理および数式

[数式→実装]にフォーカスした記事がなかったので綴ります。
今回紹介するのはCooley–Tukey型FFTのradix-2という最もシンプルなものです。

まず、DFT(離散フーリエ変換)は次の式で表されます。
$$F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) e^{-i \frac{2πtx}{N}} \qquad t=0,1,...,N-1\tag{1}$$
ここで、$N$は入力サイズで、2の累乗に限定します。
$W_N = e^{-i \frac{2π}{N}}$ とし、次のように書き換えます。
$$F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) W_{N}^{tx}\qquad t=0,1,...,N-1\tag{2}$$
入力の偶数番目の要素のみでDFTを行ったものは次の式で表されます。
$$F_{even}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k) W_\frac{N} {2}^{tk}\tag{3}$$
同様に、奇数番目の要素のみでDFTを行ったものは次の式で表されます。
$$F_{odd}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k+1) W_\frac{N}{2}^{tk}\tag{4}$$
$(3)(4)$を用いることで、$(2)$は次のように書き換えることが出来ます。
$$F(t) = F_{even}(t) +W_N^t F_{odd}(t)\tag{5}$$
ここまでは、普通にDFTを偶奇に分けて変形しただけですね。さて、ここからがFFTの重要なポイントです。実は、$F_{even}(t)，F_{odd}(t)，W_N^t$は次のような周期性を持ちます。

$F_{even}(t) = F_{even}(t+\frac{N}{2})\tag{6}$
$F_{odd}(t) = F_{odd}(t+\frac{N}{2})\tag{7}$
$W_N^t = -W_N^{t+\frac{N}{2}}\tag{8}$

$(6)(7)(8)$を$(5)$に用いることで、最終的にDFTは次のように書き換えることが出来ます。

$F(t) = F_{even}(t) +W_N^t F_{odd}(t)\qquad t=0,1,...,\frac{N}{2}-1\tag{9}$
$F(t+\frac{N}{2}) = F_{even}(t) - W_N^tF_{odd}(t)\qquad t=0,1,...,\frac{N}{2}-1\tag{10}$

$F_{even}(t)$と$F_{odd}(t)$は共通なので、これらを求めると同時に$F(t)$と$F(t+\frac{N}{2})$が求まる訳ですね。これが計算量削減の仕組みです。

以上を踏まえて実装に移りましょう。

FFTのプログラム

FFTの関数をPythonで書いてみます。
基本のアイデアはこうです。入力データを偶数要素と奇数要素に分け、$(3)(4)$のようにFFT(再帰)を行った後、$(9)(10)$を用いて0~(N/2)-1番目とN/2~N-1番目を計算します。よってFFTは再帰関数となります。

def FFT(f):
    #f:サイズNの入力データ
    N = len(f)

    if N == 1:    #Nが1のときはそのまま入力データを返す
        return f[0]

    f_even = f[0:N:2]    #fの偶数番目の要素
    f_odd = f[1:N:2]    #fの奇数番目の要素
    F_even = FFT(f_even)    #(3)偶数番目の要素でFFT
    F_odd = FFT(f_odd)    #(4)偶数番目の要素でFFT

    W_N = np.exp(-1j * (2 * np.pi * np.arange(0, N // 2)) / N)    #tが0~N/2-1番目までのWを計算した配列

    F = np.zeros(N, dtype ='complex')    #FFTの出力
    F[0:N//2] = F_even + W_N * F_odd    #(9)を計算(t:0~N/2-1)
    F[N//2:N] = F_even - W_N * F_odd    #(10)を計算(t:N/2~N-1)

    return F

動作確認

numpyのFFTと同じ出力が得られるか確認してみましょう。
次のようなパラメータで2つのsin波の合成波を作り、これを入力とします。

import numpy as np
import matplotlib.pyplot as plt

N = 64 #データ数
n = np.arange(N)
f1 = 2 #周期1
f2 = 5 #周期2
A1 = 3 #振幅1
A2 = 5 #振幅2
INPUT_f = A1 * np.sin(f1 * 2 * np.pi * (n/N)) + A2 * np.sin(f2 * 2 * np.pi * (n/N)) #2つのsin波の合成

#表示
plt.figure(figsize=(10, 5))
plt.style.use('seaborn')
plt.plot(INPUT_f)

FFTを実行します。FFTの出力は複素数のため、絶対値をかけて振幅成分を取り出し、可視化したものがこちらです。

OUTPUT_F_numpy = np.fft.fft(INPUT_f)    #numpy版
OUTPUT_F_original = FFT(INPUT_f)    #自作版

plt.figure(figsize=(10, 5))

plt.subplot(1,2,1)    #振幅成分表示(numpy版)
plt.title('numpy_ver')
plt.plot(np.abs(OUTPUT_F_numpy))

plt.subplot(1,2,2)    #振幅成分表示(自作版)
plt.title('original')
plt.plot(np.abs(OUTPUT_F_original))

無事同じ結果が得られました。

FFT(高速フーリエ変換):実装のための数式と実装例

FFTの原理および数式

FFTのプログラム

動作確認