円周率.jp - 計算量の導出

計算量の導出

ここでは FFT の計算量が $O(n\log n)$ になることについての簡単な説明を行う。

2 回の FFT $\cdots O(n\log n)$
要素ごとの積 $\cdots O(n)$
逆 FFT $\cdots O(n\log n)$
正規化 $\cdots O(n)$

という手順で行うことができるので、その計算量は $O(n\log n)$ になる。

この計算量を示す過程にはアルゴリズムによって 2 種類の導き方があるので、両者とも簡単に説明する。その際、変換元のデータは時間軸上のデータ、変換後のデータは周波数軸上のデータ、というように捉えると名前が分かりやすい。

時間間引きFFT

DFT の定義式のうち $j$ が奇数のものと偶数のものに分ける。つまり名前の通り時間軸で並べた状態から間引く形に分ける。 \[ A_k = \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} + \omega^k\sum_{j=0}^{n/2-1}a_{2j+1}\omega^{2jk} \tag{1} \]

こうなると右辺の $\sum$ の中は項数 $n/2$ の DFT になっているので、$A_k$ の方もそれに合わせる。ただしこちら側も $k$ の偶奇で分割しては意味がないので $k$ と $k+n/2$ に分けて $0\leq k \lt n/2$ にする。 \[ \begin{eqnarray} A_{k+n/2} &=& \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} + \omega^{n/2} \omega^k\sum_{j=0}^{n/2-1}a_{2j+1}\omega^{2jk}\\ &=& \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} - \omega^k\sum_{j=0}^{n/2-1}a_{2j+1}\omega^{2jk} \tag{2} \end{eqnarray} \]

ここで (1) 式と (2) 式を同時に求めることを考えると、 ${\displaystyle \sum_{j=0}^{n/2-1}} a_{2j}\omega^{2jk}$ と $\omega^k {\displaystyle \sum_{j=0}^{n/2-1}} a_{2j+1}\omega^{2jk}$ を求めればよいことが分かる。また、それらは項数 $n/2$ の DFT であるため、この方法を再帰的に適用することができることも分かる。プログラムで書くと

// x[n] の DFT を計算する時間間引き(Decimation in time) FFT
void FFT_DIT(int n, Complex* x, Complex* y) {
  if (n <= 1) return;

  Complex *x0 = x, *x1 = &x[n/2];
  Complex *y0 = y, *y1 = &y[n/2];
  for (int j = 0; j < n / 2; ++j) {
    y0[j] = x[2 * j];
    y1[j] = x[2 * j + 1];
  }
  FFT_DIT(n / 2, y0, x0);
  FFT_DIT(n / 2, y1, x1);
  for (int k = 0; k < n / 2; ++k) {
    Complex w = Complex::polar(0, -2 * M_PI * k / n);
    x0[k] = y0[k] + w * y1[k];
    x1[k] = y0[k] - w * y1[k];
  }
}

となる。なお元のデータと同サイズの一次領域 y[n] は再帰呼び出しで使いまわしている。

ここで項数 $n$ の FFT の計算量を $F(n)$ とすると、先ほどの再帰的な FFT 分割適用は

\[ F(n) = 2F(n/2) + Cn \]

となる。末尾の $Cn$ は $\omega^k$ の乗算や 2 つの FFT の結果の加減算にかかるコストを示している。この式を解くと

\[ F(n) = O(n\log n) \]

となる。

周波数間引き FFT

今度は逆に、定義式のうち周波数軸 $k$ を偶奇に間引き、時間軸 $j$ を $j$ と $j+n/2$ とに分ける。 \[ \begin{eqnarray} A_{2k} &=& \sum_{j=0}^{n/2-1} a_j\omega^{2jk} &+& \sum_{j=0}^{n/2-1} a_{j+n/2}\omega^{2jk} &=& \sum_{j=0}^{n/2-1} \left( a_j + a_{j+n/2} \right) \omega^{2jk} \tag{3}\\ A_{2k+1} &=& \sum_{j=0}^{n/2-1} a_j\omega^j\omega^{2jk} &+& \sum_{j=0}^{n/2-1} a_{j+n/2}\omega^j\omega^{n/2}\omega^{2jk} &=& \sum_{j=0}^{n/2-1} \left(a_j + \omega^{n/2}a_{j+n/2}\right)\omega^j\omega^{2jk}\\ &&&& &=& \sum_{j=0}^{n/2-1} \left( a_j - a_{j+n/2} \right) \omega^j \omega^{2jk} \tag{4}\\ \end{eqnarray} \]

ここで (3) 式と (4) 式については \[ \begin{cases} a_j^{(0)} = a_j + a_{j+n/2}\\ a_j^{(1)} = (a_j - a_{j+n/2}) \omega^j \end{cases} \]

の変換をすることでやはりサイズ $n/2$ の DFT 2 つに縮約できる。これをプログラムコードにすると

// x[n] の DFT を計算する周波数間引き(Decimation in frequency) FFT
void FFT_DIF(int n, Complex* x, Complex* y) {
  if (n <= 1)
    return;

  Complex *y0 = y, *y1 = &y[n/2];
  for (int j = 0; j < n / 2; ++j) {
    Complex w = Complex::polar(0, -2 * M_PI * j / n);
    y0[j] = x0[j] + x1[j];
    y1[j] = (x0[j] - x1[j]) * w;
  }
  FFT_DIF(n / 2, y0, x0);
  FFT_DIF(n / 2, y1, x1);
  for (int k = 0; k < n / 2; ++k) {
    x[2*k  ] = y0[k];
    x[2*k+1] = y1[k];
  }
}

となる。