FFTを数式で書いてみた - 情報系の手考ノート

FFTについてgoogle先生にお聞きしてみると，実装できるところまで解説したらおしまい，という形をとってるものが多い気がします。実際，実装できればいいというのはそのとおりだと思いますが，やはりアルゴリズムの理解という意味では自分で手を動かして計算してみるのが一番ではないかと思います。

というわけでFFTに対する理解を深めるために，式を立ててみました。自分へのメモがてらまとめておきます。

ちなみにデータ数 $N$ が2の累乗数のときに限定して話を進めます。

なお，自分で考えて導出した式なので間違ってる可能性も大いにあります。

※2019/5/9 : 順変換の式の $x$ の範囲が一部間違ってたので修正しました。 ※2019/5/10 : 分割した式の $w$ の添字が一部間違ってたので修正しました。

FFTとは

FFTとは Fast Fourier Transform 略で，離散フーリエ変換を計算機上で高速に計算するアルゴリズムです。離散フーリエ変換の計算量は $O(N^ 2)$ ですが，FFTの計算量はデータ数 $N$ が2の累乗数のときに限り $O(N \log N)$ です。

離散フーリエ変換(Descrete Fourier Transform)

離散フーリエ変換(以下，DFT)はフーリエ変換を離散化したもので，周波数解析などでよく使われます。順変換は，虚数単位 $i$ ，データ数 $N$ ，ネイピア数 $e$ ，円周率 $\pi$ を用いて

$\begin{aligned} F(x) = \sum_{t=0}^{N-1} f(t) e^{-i \frac{2 \pi tx}{N}} \end{aligned}$

と定義され，逆変換は

$\begin{aligned} f(t) = \frac{1}{N} \sum_{x=0}^{N-1} F(x) e^{i \frac{2 \pi xt}{N}} \end{aligned}$

と定義されます。

この式を変形していくと，データ数 $N$ のDFTはデータ数 $\frac{N}{2}$ のDFT2に分割できます。この分割作業を再帰的に繰り返していくことで計算量を $O(N \log N)$ にまで落とすことができるのです。

Cooley-Tukey型アルゴリズム

今回は，いくつかあるらしいFFTのアルゴリズムのうち，Cooley-Tukey型と言われるアルゴリズムに従って式を変形していきます。

まず，計算を簡単化するためにデータ数は自然数 $k$ を用いて $N = 2^ k$ とかけるとします。それと，DFTの式を

$\begin{aligned} F(x) = \sum_{t=0}^{N-1} f(t) w^{tx}_{N} \end{aligned}$

$\begin{aligned} f(t) = \frac{1}{N} \sum_{x=0}^{N-1} F(x) w^{-xt}_{N} \end{aligned}$

と書き換えます。このとき， $e^{-i \frac{2 \pi tx}{N}} = w^{tx}_{N}$ としました。

順変換

まず，順変換について考えます。順変換の式を変形すると

$\begin{aligned} F(x) &= \sum_{t=0}^{N-1} f(t) w^{tx}_{N} \\ &= \sum_{t=0}^{\frac{N}{2}-1} [ f(2t) w^{2t \cdot x}_{N} + f(2t + 1) w^{(2t + 1) \cdot x}_{N} ] \\ &= \sum_{t=0}^{\frac{N}{2}-1} f(2t) w^{tx}_{\frac{N}{2}} + w^{x}_{N} \sum_{x=0}^{\frac{N}{2}-1} f(2t + 1) w^{tx}_{\frac{N}{2}} \end{aligned}$

と書くことができます。この式は $x$ の範囲を $0 \le x \le \frac{N}{2} - 1$ とすると

$\begin{aligned} F(x) = \sum_{t=0}^{\frac{N}{2}-1} f(2t) w^{tx}_{\frac{N}{2}} + w^{x}_{N} \sum_{t=0}^{\frac{N}{2}-1} f(2t + 1) w^{tx}_{\frac{N}{2}} \end{aligned}$

$\begin{aligned} F(x + \frac{N}{2}) = \sum_{t=0}^{\frac{N}{2}-1} f(2t) w^{t\left(x + \frac{N}{2}\right)}_{\frac{N}{2}} + w^{x + \frac{N}{2}}_{N} \sum_{t=0}^{\frac{N}{2}-1} f(2t + 1) w^{t\left(x + \frac{N}{2}\right)}_{\frac{N}{2}} \end{aligned}$

という2つの式であらわすことができます。ここで， $w$ は指数関数の性質を考えると

$\begin{aligned} w^{x + \frac{N}{2}}_{N} = -w^{x}_{N} \end{aligned}$

$\begin{aligned} w^{x + N}_{N} = w^{x}_{N} \end{aligned}$

という性質を持ちます。これを考えると，上の2つの式はそれぞれ

$\begin{aligned} F(x) = \sum_{t=0}^{\frac{N}{2}-1} f(2t) w^{tx}_{\frac{N}{2}} + w^{x}_{N} \sum_{t=0}^{\frac{N}{2}-1} f(2t + 1) w^{tx}_{\frac{N}{2}} \end{aligned}$

$\begin{aligned} F(x + \frac{N}{2}) = \sum_{t=0}^{\frac{N}{2}-1} f(2t) w^{tx}_{\frac{N}{2}} - w^{x}_{N} \sum_{t=0}^{\frac{N}{2}-1} f(2t + 1) w^{tx}_{\frac{N}{2}} \end{aligned}$

として書くことができます。この式を見てみると，データ数 $\frac{N}{2}$ のDFTが2つあらわれていることがわかります。

では，2つのDFTを

$\begin{aligned} f_{\frac{N}{2}}(x) = \sum_{t=0}^{\frac{N}{2}-1} f(2t) w^{tx}_{\frac{N}{2}} \end{aligned}$

$\begin{aligned} f_{\frac{N}{2}}(x + \frac{N}{2}) = \sum_{t=0}^{\frac{N}{2}-1} f(2t + 1) w^{tx}_{\frac{N}{2}} \end{aligned}$

として，これら2つを分割することを考えてみます。最初に分割したときと同様の手順を取ることで分割してみると，

$\begin{aligned} f_{\frac{N}{2}}(x + \frac{N}{2}b_0) &= \sum_{t=0}^{\frac{N}{4}-1} f(4t + 2b_0) w^{tx}_{\frac{N}{4}} + w^x_{\frac{N}{2}} \sum_{t=0}^{\frac{N}{4}-1} f(4t + 2b_0 + 1) w^{tx}_{\frac{N}{4}} \end{aligned}$

$\begin{aligned} f_{\frac{N}{2}}(x + \frac{N}{2}b_0 + \frac{N}{4}) &= \sum_{t=0}^{\frac{N}{4}-1} f(4t + 2b_0) w^{tx}_{\frac{N}{4}} - w^x_{\frac{N}{2}} \sum_{t=0}^{\frac{N}{4}-1} f(4t + 2b_0 + 1) w^{tx}_{\frac{N}{4}} \end{aligned}$

となります。ここで， $x$ は $0 \le x \le \frac{N}{4} - 1$ ， $b_0 \in \{0, 1 \}$ です。

このような手順を繰り返していくと，最終的にデータ数 $N$ のDFTの順変換は1以上 $k$ 以下の整数 $l$ を用いて，

$\begin{aligned} F(x + \frac{N}{2} b_0) = f_{\frac{N}{2}}(x) + (-1)^{b_0} w^{x}_{N} f_{\frac{N}{2}}(x + \frac{N}{2})) \end{aligned}$

$\begin{aligned} f_{\frac{N}{2^ l}}(x + \frac{N}{2^{l+1}}(2^ l b_0 + \cdots + 2b_{l-1})) = \sum^{\frac{N}{2^{l+1}}-1}_{t=0} f(2^ l t + 2^{l-1} b_{l-1} + \cdots + b_0) w^{tx}_{\frac{N}{2^ l}} \\ (0 \le x \le \frac{N}{2^ l}-1) \end{aligned}$

と置くと $l$ が $1 \le l \le k-1$ であり，かつ $0 \le x \le \frac{N}{2^{l+1}}$ としたとき $f_{\frac{N}{2^ l}}$ は

$\begin{aligned} f_{\frac{N}{2^ l}}(x + \frac{N}{2^{l+1}}(2^ l &b_0 + \cdots + 2b_{l-1} + b_l)) \\ =& f_{\frac{N}{2^{l+1}}}(x + \frac{N}{2^{l+1}}(2^ l b_0 + \cdots + 2b_{l-1})) \\ &+ (-1)^{b_l} w^{x}_{\frac{N}{2^{l}}} f_{\frac{N}{2^{l+1}}}(x + \frac{N}{2^{l+1}}(2^ l b_0 + \cdots + 2b_{l-1} + 1)) \end{aligned}$

と書くことができます。

これで順変換の式を求めることができました。

逆変換

次に，逆変換について考えます。逆変換の式中にある $w^{-xt}_{N}$ を $\omega^{xt}_{N}$ と置き換えてみると

$\begin{aligned} f(t) = \frac{1}{N} \sum_{x=0}^{N-1} F(x) \omega^{xt}_{N} \end{aligned}$

と書き直すことができます。そして $\omega$ は $w$ と同様に

$\begin{aligned} \omega^{x + \frac{N}{2}}_{N} = -w^{x}_{N} \end{aligned}$

$\begin{aligned} \omega^{x + N}_{N} = w^{x}_{N} \end{aligned}$

という性質を持ちます。つまり，順変換のときとまったく同じように式変形することができます。順変換と同脳の変形を行うと逆変換は

$\begin{aligned} f(t + \frac{N}{2} b_0) = \frac{1}{N} (F_{\frac{N}{2}}(t) + (-1)^{b_0} w^{-t}_{N} F_{\frac{N}{2}}(t + \frac{N}{2})) \end{aligned}$

$\begin{aligned} F_{\frac{N}{2^ l}}(t + \frac{N}{2^{l+1}}(2^ l b_0 + \cdots + 2b_{l-1})) = \sum^{\frac{N}{2^{l+1}}-1}_{x=0} F(2^ l x + 2^{l-1} b_{l-1} + \cdots + b_0) w^{-xt}_{\frac{N}{2^ l}} \\ (0 \le t \le \frac{N}{2^ l}-1) \end{aligned}$

と置くと $l$ が $1 \le l \le k-1$ であり，かつ $0 \le t \le \frac{N}{2^{l+1}}$ としたとき

$\begin{aligned} F_{\frac{N}{2^ l}}(t + \frac{N}{2^{l+1}}(2^ l &b_0 + \cdots + 2b_{l-1} + b_l)) \\ =& F_{\frac{N}{2^{l+1}}}(t + \frac{N}{2^{l+1}}(2^ l b_0 + \cdots + 2b_{l-1})) \\ &+ (-1)^{b_l} w^{-t}_{\frac{N}{2^{l}}} F_{\frac{N}{2^{l+1}}}(t + \frac{N}{2^{l+1}}(2^ l b_0 + \cdots + 2b_{l-1} + 1)) \end{aligned}$

となります。

総括

まとめです。

DFT前の関数を $f$ ，DFT後の関数を $F$ としたとき，データ数 $N = 2^ k$ のDFTの順変換は $0 \le x \le \frac{N}{2}-1$ として

$\begin{aligned} F(x + \frac{N}{2} b_0) = f_{\frac{N}{2}}(x) + (-1)^{b_0} w^{x}_{N} f_{\frac{N}{2}}(x + \frac{N}{2}) \end{aligned}$

と書くことができ， $1 \le l \le k-1$ であるような $l$ を用いて $f_{\frac{N}{2^ l}}$ は $0 \le x \le \frac{N}{2^{l+1}}-1$ の範囲で

$\begin{aligned} f_{\frac{N}{2^ k}}(2^{k-1} b_0 + \cdots + b_{k-1}) = f(2^{k-1} b_{k-1} + \cdots + b_0) \end{aligned}$

と書くことができます。

同様にDFTの逆変換は

$\begin{aligned} f(t + \frac{N}{2} b_0) = \frac{1}{N} (F_{\frac{N}{2}}(t) + (-1)^{b_0} w^{-t}_{N} F_{\frac{N}{2}}(t + \frac{N}{2})) \end{aligned}$

と書くことができ， $F_{\frac{N}{2^ l}}$ は

$\begin{aligned} F_{\frac{N}{2^ k}}(2^{k-1} b_0 + \cdots + b_{k-1}) = F(2^{k-1} b_{k-1} + \cdots + b_0) \end{aligned}$

と書くことができます。

DFTの順変換，逆変換の式中で用いた $b_0, b_1, \cdots , b_{k-1}$ は1か0をとるものとします。

上の式の $f_{\frac{N}{2^ k}}$ ， $F_{\frac{N}{2^ k}}$ を見てみると，それぞれビットの順序を入れ返す作業が行われてることがわかります。よくFFTの解説で見るビットの順序を入れ替える操作はこの式と対応づいているわけです。

実際に実装する場合，0から $k-1$ まで $f_{\frac{N}{2^ l}}$ を計算することになります。これは $k$ 回の計算なので $\log N$ のオーダーです。さらに各 $f_{\frac{N}{2^ l}}$ は $N$ のオーダーで計算できるので最終的に，計算量が $O(N \log N)$ になります。

ちなみに，多倍長演算にもFFTと同じ考え方が使われているらしいので，多倍長演算用のアルゴリズムなんかも調べてみたいです。