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

##概要
スペクトル法に利用するFFTについて，self-soating型[^ishioka2004][^yamamoto2003]のアルゴリズムを調べた．Fortran90でコーディングし，直接計算（DFT）と比較してデバッグした．[](最後に参考文献[^ishioka2004]の著者が公開している数値計算ライブラリISPACK[^ispack]のFFTと実行速度の比較を行った．)

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

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

##本記事の目的

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

##本記事の目的でないもの

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

[^fftw]: FFTW (http://www.fftw.org/) 等

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

##定義

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

* **DFT（逆変換）**

```math
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（正変換）**

```math
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型[^yamamoto2003]
　入力配列=出力配列(in place)で省メモリ．Wikipedia[^wikipedia]等，解説しているサイトはたくさんあるので省略．出力を入力順にするにはビット反転操作が必要．

[^wikipedia]: Wikipedia : 高速フーリエ変換 (https://ja.wikipedia.org/wiki/高速フーリエ変換)

* Stockham型[^yamamoto2003]
　入力配列の他に同サイズの作業用配列が必要(out of place)．出力は入力順に対応する．self-soating型とも呼ばれる．文献[^ishioka2004]に解説があるため，今回はStockham型アルゴリズムを確認する．[](データ長は高々$2^{10}$程度を想定するため省メモリ化をさほど気にする必要がなく（むしろデータ長はキャッシュに乗り切るかどうかが大切なはずだが浅学のため検証していない），)

###Stockham型の複素FFTアルゴリズム
　式(1)を行列形式で書くと

```math
\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}$に対して，

```math
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$のとき

```math
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$のとき

```math
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))の左辺は，

```math
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$を用いて表すと

```math
\hat{s}=R\hat{q}+\hat{r}\\
s=Qr+q
```
と書ける．よって式(4)で定義された$F_S$は

```math
\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}
```

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

```math
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$
>&nbsp;&nbsp;&nbsp;&nbsp;$F_R(\hat{p})=f(\hat{p})$
>end do
>$P=N$
>do $j=1,l$
>&nbsp;&nbsp;&nbsp;&nbsp;$P=P/Q_j$
>&nbsp;&nbsp;&nbsp;&nbsp;$Q=Q_j$
>&nbsp;&nbsp;&nbsp;&nbsp;$R=N/(PQ)$
>&nbsp;&nbsp;&nbsp;&nbsp;do $\hat{r}=0,R-1$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;do $\hat{p}=0,P-1$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$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})$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;end do
>&nbsp;&nbsp;&nbsp;&nbsp;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)および疑似コードは文献[^ishioka2004]から一部改変して引用した．

#### 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)$
>&nbsp;&nbsp;&nbsp;&nbsp;$P=P/Q_1=8/2=4$
>&nbsp;&nbsp;&nbsp;&nbsp;$Q=Q_1=2$
>&nbsp;&nbsp;&nbsp;&nbsp;$R=N/(PQ)=8/(4\times 2)=1$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(0)=F_1(0)+F_1(4)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(1)=F_1(0)-F_1(4)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(2)=F_1(1)+F_1(5)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(3)=\left\\{F_1(1)-F_1(5)\right\\}\exp\left(\frac{\pi}{4}i\right)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(4)=F_1(2)+F_1(6)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(5)=\left\\{F_1(2)-F_1(6)\right\\}\exp\left(\frac{\pi}{2} i\right)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(6)=F_1(3)+F_1(7)$
>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{2\;1}(7)=\left\\{F_1(3)-F_1(7)\right\\}\exp\left(\frac{3\pi}{4}i\right)$

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

```math
\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$は次のようになる．

```math
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を文献[^ishioka2004]を参照してまとめた．Stockham型のアルゴリズムを行列形式に書き下すと，密な係数行列を疎行列に分解する事で効率的な演算を達成することがわかった．本記事では触れなかったが，Cooley-Tukey型のバタフライ演算も同様に考えられる[^wikipedia]．
　Fortran90でコーディングしたプログラムについては，次回以降で性能評価を行う予定である．[](自作のコードを公開されているライブラリ[^ispack]と比較した結果，実行時間は数倍のオーダーで収まった．)当然のことだが，速度と正しさ（十分検証されていること）を求めるならFFTW[^fftw]，Intel MKL[^mkl]，ISPACK[^ispack]等を使用するべきである．

[^ispack]: ISPACK: 科学計算のためのFORTRANライブラリ (https://www.gfd-dennou.org/library/ispack/)
[^mkl]: Intel Math Kernel Library (https://software.intel.com/en-us/mkl)

##次回以降の予定
1. [実FFTのアルゴリズム](https://qiita.com/toya42/items/dacf940f16334e0d8543)
2. ライブラリ（Intel MKL, ISPACK）の利用方法（+自作コードの性能評価）