趣味でスペクトル法_#001_FFTの基礎 #Fortran

概要

スペクトル法に利用するFFTについて，self-soating型¹²のアルゴリズムを調べた．Fortran90でコーディングし，直接計算（DFT）と比較してデバッグした．

更新履歴

2019.06.16 : 公開．
2019.06.16 : 小さなミスを修正
2019.06.22 : 次の記事のリンクを追加．"定義"に記号に関する注意書きを追加

本記事の目的

FFTのアルゴリズムを理解し，Fortran90でコーディングする．
QiitaおよびGithubの使い方の確認．

本記事の目的でないもの

高速あるいは省メモリなFFTライブラリ³の作成
数学的に厳密なDFTの解説
複数のFFTアルゴリズムの比較検討

略記

DFT: 離散フーリエ変換 (discrete Fourier Transform)
FFT: 高速フーリエ変換 (fast Fourier Transform)

定義

通常，波数空間で定義される値に$F,k$を用い，物理空間で定義される値に$f,x$を用いるが，本記事では逆になっている．（特に理由はないが支障もない）

DFT（逆変換）

F\left(\hat{k}\right)=\sum_{x=0}^{N-1} f\left(x\right) \exp\left(2\pi i \frac{x\cdot\hat{k}}{N}\right) \tag{1}\\
%\left(k=0,1,...,N-1\right)

DFT（正変換）

f\left(x\right)=\frac{1}{N}\sum_{\hat{k}=0}^{N-1} F\left(\hat{k}\right) \exp\left(-2\pi i \frac{\hat{k}\cdot x}{N}\right) \tag{2}\\
=\frac{1}{N}\overline{ \sum_{\hat{k}=0}^{N-1} \overline{F\left(\hat{k}\right)} \exp\left(2\pi i \frac{x\cdot\hat{k}}{N}\right)}\\
%\left(x=0,1,...,N-1\right)\\

入力データ長と出力データ長は等しい．$\left(x,k=0,1,...,N-1\right)$．$f,F$は複素数データ列である．また複素数$C$の複素共役は，$C^*$あるいは$\bar{C}$と書く．

FFTアルゴリズム

上で示したように正変換は逆変換を利用して得られるため，逆変換のみを考えれば良い．

Cooley-Tukey型²
　入力配列=出力配列(in place)で省メモリ．Wikipedia⁴等，解説しているサイトはたくさんあるので省略．出力を入力順にするにはビット反転操作が必要．
Stockham型²
　入力配列の他に同サイズの作業用配列が必要(out of place)．出力は入力順に対応する．self-soating型とも呼ばれる．文献¹に解説があるため，今回はStockham型アルゴリズムを確認する．

Stockham型の複素FFTアルゴリズム

　式(1)を行列形式で書くと

\left(\begin{matrix}
F(0)  \\
\vdots\\
F\left(\hat{k}\right)\\
\vdots\\
F(N-1)
\end{matrix}\right)
=
\left(\begin{matrix}
(\hat{k},x)成分が\\\\
\exp\left(2\pi i \frac{x\cdot\hat{k}}{N}\right)なる\;\;\\\\
N\times N行列
%g(0,0) & \ldots & g(0,x) & \ldots & g(0,N-1) \\
%\vdots &  \ddots & \vdots& &\vdots\\
%g(\hat{k},0) &\ldots &g(\hat{k},x)& \ldots & g(\hat{k},N-1) \\
%\vdots & \\
%g(N-1,0) & \ldots  & g(N-1,x) &\ldots & g(N-1,N-1)
\end{matrix}\right)
\left(\begin{matrix}
f(0)  \\
\vdots\\
f(x) \\
\vdots\\
f(N-1)
\end{matrix}\right) \tag{3}

であり，DFTは$N\times N$の正方行列$A$と長さ$N$の列ベクトルの積なので，その計算量は$\mathcal{O}(N^2)$となる．

漸化式の導入

　さてここで，$N=PS,\;0\leq\hat{p}\leq P-1,\;0\leq\hat{s}\leq S-1$なる正の整数$P,S$，及び整数$\hat{p},\hat{s}$に対して，

F_{s}\left(S\hat{p}+\hat{s}\right)
=\exp\left(2\pi i \frac{\hat{p}\hat{s}}{PS}\right)
\sum_{s=0}^{S-1}f\left(Ps+\hat{p}\right)
\exp\left(2\pi i \frac{\hat{s}s}{S}\right) \tag{4}

なる$X_s$を導入する．すると，

$S=1,P=N$のとき

F_1\left(\hat{p}\right)=f\left(\hat{p}\right),\;\;
\left(\hat{p}=0,1,...,P-1\right) \tag{5}

$S=N,P=1$のとき

F_N\left(\hat{s}\right)=
\sum_{s=0}^{S-1}f\left(s\right)
\exp\left(2\pi i \frac{\hat{s}s}{S}\right)
=F\left(\hat{s}\right)
,\;\;
\left(\hat{s}=0,1,...,S-1\right) \tag{6}

となる．つまりFFTで求めたい式(1)（あるいは(3))の左辺は，

F\left(\hat{k}\right)=F_N\left(\hat{k}\right)

で得られることになる．
式(5)より$F_1$は既知なので，$F_S$の$S$を1から$N$まで増す漸化式を作ればよい．

　$N=PS$に加えて$S=QR$とおき($P,Q,R,S$は正の整数），$0\leq s,\hat{s}\leq S-1$なる $s,\hat{s}$ を$0\leq q,\hat{q}\leq Q-1$と$0\leq r,\hat{r}\leq R-1$を用いて表すと

\hat{s}=R\hat{q}+\hat{r}\\
s=Qr+q

と書ける．よって式(4)で定義された$F_S$は

\begin{align}
F_S\left(S\hat{p}+\hat{s}\right)&=F_{QR}\left(QR\hat{p}+R\hat{q}+\hat{r}\right)\\
&=
\exp\left(2 \pi i \frac{\hat{p}\left(R\hat{q}+\hat{r}\right)}{PQR}\right)
\sum_{q=0}^{Q-1}\sum_{r=0}^{R-1}f\left(P\left(Qr+q\right)+\hat{p}\right)
\exp\left(2\pi i \left\{\frac{q\hat{q}}{Q}+\frac{r\hat{r}}{R}\right\}\frac{q\hat{r}}{QR}\right)\\
&=
\exp\left(2 \pi i \frac{\hat{p}\hat{q}}{PQ}\right)
\sum_{q=0}^{Q-1}
\exp\left(2 \pi i \frac{\hat{q}q}{Q}\right)
\exp\left(\frac{\left(Pq+\hat{p}\right)\hat{r}}{PQR}\right)
\sum_{r=0}^{R-1}
f\left(PQr+\left(Pq+\hat{p}\right)\right)
\exp\left(2\pi i \frac{r\hat{r}}{R}\right)\\
&=
\exp\left(2 \pi i \frac{\hat{p}\hat{q}}{PQ}\right)
\sum_{q=0}^{Q-1}
\exp\left(2 \pi i \frac{\hat{q}q}{Q}\right)
F_{R}\left(R\left(Pq+\hat{p}\right)+\hat{r}\right) \tag{7}
\end{align}

よって漸化式は次のようにまとめられる．

F_{QR}\left(QR\hat{p}+R\hat{q}+\hat{r}\right)
=
\exp\left(2 \pi i \frac{\hat{p}\hat{q}}{PQ}\right)
\sum_{q=0}^{Q-1}
\exp\left(2 \pi i \frac{\hat{q}q}{Q}\right)
F_{R}\left(R\left(Pq+\hat{p}\right)+\hat{r}\right) \tag{8}

疑似コード

　疑似コードは，$N=\prod_{j=1}^{l}Q_j$とすると，次のように書ける．

do $\hat{p}=0,N-1$
    $F_R(\hat{p})=f(\hat{p})$
end do
$P=N$
do $j=1,l$
    $P=P/Q_j$
    $Q=Q_j$
    $R=N/(PQ)$
    do $\hat{r}=0,R-1$
        do $\hat{p}=0,P-1$
            $F_{QR}(QRp+R\hat{q}+\hat{r})=\exp\left(2\pi i \frac{\hat{p}\hat{q}}{PQ}\right)\sum_{q=0}^{Q-1}\exp\left(2\pi i \frac{\hat{q}q}{Q}\right)F_R(PRq+R\hat{p}+\hat{r})$
        end do
    end do
end do

$F_R\rightarrow F_{QR}$の変換に必要な演算量は，疑似コードから$RPQ^2=NQ$である．
よって，$F_1\rightarrow F_{Q_1}\rightarrow F_{Q_1 Q_2} \rightarrow ...\rightarrow F_N=F$の変換に必要な演算量は$N\left(Q_1+Q_2+...Q_l\right)$と計算できる．特に$N=2^l$のときは，$N(2\times l) =2N\log_2N$となる．密行列の積を直接計算するDFTの演算量が$\mathcal{O}(N^2)$だったことと比較すると，FFTにより非常に効率よく変換できるようになった．

なお，式(3)-(8)および疑似コードは文献¹から一部改変して引用した．

N=8の場合

　このままではわかりづらいため，小さい$N$の場合について書き下してみる．

$N=8\;\left(Q_1=Q_2=Q_3=2\right)$のとき

$F_1\rightarrow F_2$の操作を疑似コードで表すと，

$F_1(0:7)=f(0:7)$
$P=8$
$\left(j=1\right)$
    $P=P/Q_1=8/2=4$
    $Q=Q_1=2$
    $R=N/(PQ)=8/(4\times 2)=1$
            $F_{2\;1}(0)=F_1(0)+F_1(4)$
            $F_{2\;1}(1)=F_1(0)-F_1(4)$
            $F_{2\;1}(2)=F_1(1)+F_1(5)$
            $F_{2\;1}(3)=\left\{F_1(1)-F_1(5)\right\}\exp\left(\frac{\pi}{4}i\right)$
            $F_{2\;1}(4)=F_1(2)+F_1(6)$
            $F_{2\;1}(5)=\left\{F_1(2)-F_1(6)\right\}\exp\left(\frac{\pi}{2} i\right)$
            $F_{2\;1}(6)=F_1(3)+F_1(7)$
            $F_{2\;1}(7)=\left\{F_1(3)-F_1(7)\right\}\exp\left(\frac{3\pi}{4}i\right)$

となる．上の演算を行列形式で書けば

\left(\begin{matrix}
F_2(0)\\
F_2(1)\\
F_2(3)\\
F_2(4)\\
F_2(5)\\
F_2(6)\\
F_2(7)
\end{matrix}\right)
=
\left(\begin{array}{ccccccc}
1&0&0&0&1&0&0&0\\
1&0&0&0&-1&0&0&0\\
0&1&0&0&0&1&0&0\\
0&e^{\frac{\pi}{4} i}&0&0&0&-e^{\frac{\pi}{4} i}&0&0\\
0&0&1&0&0&0&1&0\\
0&0&i&0&0&0&-i&0\\
0&0&0&1&0&0&0&1\\
0&0&0&e^{\frac{3\pi}{4} i}&0&0&0&-e^{\frac{3\pi}{4} i}
\end{array}\right)
\left(\begin{matrix}
F_1(0)\\
F_1(1)\\
F_1(3)\\
F_1(4)\\
F_1(5)\\
F_1(6)\\
F_1(7)
\end{matrix}\right)

非零要素が$2N\left(=16 \right)$個の疎行列（$A_2$とする）と長さ$N$の列ベクトルの掛け算となり，演算量は確かに$2N$であることがわかる．$F_4=F_{2\;2}\leftarrow F_2$の変換の係数行列$A_4$および，$F_8=F_{2\;4}\leftarrow F_4$の変換の係数行列$A_8$は次のようになる．

A_4=
\left(\begin{array}{ccccccc}
1&0&0&0&1&0&0&0\\
0&1&0&0&0&1&0&0\\
1&0&0&0&-1&0&0&0\\
0&1&0&0&0&-1&0&0\\
0&0&1&0&0&0&1&0\\
0&0&0&1&0&0&0&1\\
0&0&1&0&0&0&-1&0\\
0&0&0&1&0&0&0&1
\end{array}\right)\\
\\
A_8=
\left(\begin{array}{ccccccc}
1&0&0&0&1&0&0&0\\
0&1&0&0&0&1&0&0\\
0&0&1&0&0&0&1&0\\
0&0&0&1&0&0&0&1\\
1&0&0&0&-1&0&0&0\\
0&1&0&0&0&-1&0&0\\
0&0&1&0&0&0&-1&0\\
0&0&0&1&0&0&0&-1
\end{array}\right)

ここで確かに$A=A_8\; A_4\; A_2$である（maximaで検算した）．
　つまり，FFTは$N\times N$の密行列$A$を各々非零要素が$2N$個の疎行列に分解して右から順に（零要素を無視して）計算していくことで，効率的に$F_1\rightarrow F_N$を達成するアルゴリズムと言える．$N=2^l$のとき$A$は$l$個の行列に分解できるので，演算量は$2N\times l=2N\log_2N$と見積もることができる．

まとめ

　勉強のため，DFTの高速計算アルゴリズムであるFFTを文献¹を参照してまとめた．Stockham型のアルゴリズムを行列形式に書き下すと，密な係数行列を疎行列に分解する事で効率的な演算を達成することがわかった．本記事では触れなかったが，Cooley-Tukey型のバタフライ演算も同様に考えられる⁴．
　Fortran90でコーディングしたプログラムについては，次回以降で性能評価を行う予定である．当然のことだが，速度と正しさ（十分検証されていること）を求めるならFFTW³，Intel MKL⁵，ISPACK⁶等を使用するべきである．

次回以降の予定

実FFTのアルゴリズム
ライブラリ（Intel MKL, ISPACK）の利用方法（+自作コードの性能評価）

石岡圭一(2004)『スペクトル法による数値計算入門』東京大学出版 ↩
電通大山本有作教授の講義資料(2003)『高速フーリエ変換とその並列化(I)』(http://na.scitec.kobe-u.ac.jp/~yamamoto/lectures/parallelFFT/parallelFFT1.PDF) 2019.6.9閲覧 ↩
FFTW (http://www.fftw.org/) 等 ↩
Wikipedia : 高速フーリエ変換 (https://ja.wikipedia.org/wiki/高速フーリエ変換) ↩
Intel Math Kernel Library (https://software.intel.com/en-us/mkl) ↩
ISPACK: 科学計算のためのFORTRANライブラリ (https://www.gfd-dennou.org/library/ispack/) ↩