FFTWで1次元実数データをフーリエ変換

[**FFTW**](http://www.fftw.org/)は、[高速フーリエ変換](https://ja.wikipedia.org/wiki/高速フーリエ変換)（FFT）に基づいて[離散フーリエ変換](https://ja.wikipedia.org/wiki/離散フーリエ変換)（DFT）を高速に計算できるライブラリのひとつ。任意の長さのデータが扱え、多次元FFTや実数に最適化した変換も用意されている。

これまでスペクトル解析や微分など様々な計算でお世話になっていたが、毎回使い方を忘れて定数倍ずれたりしていたので復習することにした。

## 実行環境など
記事では以下の環境を用いている。

- FFTW 3.3.8
- 言語: Fortran90 ← ["Calling FFTW from Legacy Fortran"](http://www.fftw.org/fftw3_doc/Calling-FFTW-from-Legacy-Fortran.html)
- コンパイラ: gcc (gfortran含む)
- OS: Ubuntu 16.04

サンプルコードは以下の点を考慮している。

- データの長さ `N` が偶数・奇数のどちらでも適切に動く ※ `N/2` は小数点以下切り捨て
- 配列の添字は必ず0から始める

## TL;DR: プラン作成一覧
FFTの種類に応じて、サブルーチン名や引数、配列の型と大きさがどうなっていればいいかを纏める。

```fortran:plan_list.f90
program plan_list
  implicit none
  include 'fftw3.f'

  integer, parameter :: N = 1024
  integer, parameter :: FLAGS = ior(FFTW_MEASURE, FFTW_DESTROY_INPUT)

  integer(8) :: plan
  real(8)    :: r_time(0:N-1), r_freq(0:N-1)
  complex(8) :: c_time(0:N-1), c_freq(0:N-1), cp_freq(0:N/2)

  ! complex <-> complex
  call dfftw_plan_dft_1d(plan, N, c_time, c_freq, FFTW_FORWARD , FLAGS)
  call dfftw_plan_dft_1d(plan, N, c_freq, c_time, FFTW_BACKWARD, FLAGS)
  !! (execute)
  call dfftw_execute_dft(plan, c_time, c_freq)

  ! real <-> complex(half)
  call dfftw_plan_dft_r2c_1d(plan, N, r_time , cp_freq, FLAGS)
  call dfftw_plan_dft_c2r_1d(plan, N, cp_freq, r_time , FLAGS)
  !! (execute)
  call dfftw_execute_dft_r2c(plan, r_time , cp_freq)
  call dfftw_execute_dft_c2r(plan, cp_freq, r_time )

  ! real <-> real
  !! real <-> half complex
  call dfftw_plan_r2r_1d(plan, N      , r_time, r_freq, FFTW_R2HC   , FLAGS)
  call dfftw_plan_r2r_1d(plan, N      , r_freq, r_time, FFTW_HC2R   , FLAGS)
  !! DCT I-IV and DST I-IV
  call dfftw_plan_r2r_1d(plan, 2*(N-1), r_time, r_freq, FFTW_REDFT00, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*N    , r_time, r_freq, FFTW_REDFT10, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*N    , r_time, r_freq, FFTW_REDFT01, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*N    , r_time, r_freq, FFTW_REDFT11, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*(N+1), r_time, r_freq, FFTW_RODFT00, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*N    , r_time, r_freq, FFTW_RODFT10, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*N    , r_time, r_freq, FFTW_RODFT01, FLAGS)
  call dfftw_plan_r2r_1d(plan, 2*N    , r_time, r_freq, FFTW_RODFT11, FLAGS)
  !! DHT
  call dfftw_plan_r2r_1d(plan, N      , r_time, r_freq, FFTW_DHT    , FLAGS)
  !! (execute)
  call dfftw_execute_r2r(plan, r_time, r_freq)

  ! (destroy)
  call dfftw_destroy_plan(plan)

  stop
end program
```

<!--
単精度の場合は、サブルーチンは `dfftw_...` でなく `sfftw_...` を使う。
-->

## FFTWの使用準備
既にコンパイルまでできるなら飛ばしていい。

### インストール
ここではホームディレクトリ以下にインストールする。

- 管理者権限が必要なところへインストールする場合は、最後の `make install` だけ管理者権限で実行すればいい
- コンパイラやコンパイルオプションなどを指定したい場合は、configureのところで指定する

```shell-session:install_FFTW
$ wget http://www.fftw.org/fftw-3.3.8.tar.gz
$ tar zxf fftw-3.3.8.tar.gz
$ cd fftw-3.3.8/
$ ./configure --prefix=${HOME}/local
$ make
$ make install
```

```shell-session:check_files
$ ls ${HOME}/local/include  #=> fftw3.f
$ ls ${HOME}/local/lib      #=> libfftw3.a
```

### コンパイルとリンク

- プログラムにはFortran用のincludeファイル `fftw3.f` を読み込ませる
	- 拡張子からわかる通りFORTRAN77形式だが、Fortran90以降でも問題なく使える
- コンパイル時にincludeファイルの場所を教える
- リンク時にライブラリ `libfftw3.a` の場所を教え、**使用を明示**する
	- コンパイルとリンクを同時にする場合は両方の指定が必要
	- 環境変数の設定によっては `-I`, `-L` は不要

```fortran:sample.f90
program sample
  implicit none
  include 'fftw3.f'

  write(*, *) FFTW_ESTIMATE

  stop
end program sample
```

```shell-session:compile_and_link
$ f95 -c sample.f90 -I${HOME}/local/include
$ f95    sample.o   -L${HOME}/local/lib -lfftw3
$ ls
a.out  sample.f90  sample.o
$ ./a.out
          64
```

## 複素数データのFFT
最も基本的な使い方。効率を気にしないならこれだけでも行ける。

```fortran:fftw_c2c.f90
program fftw_c2c
  implicit none
  include 'fftw3.f'

  real(8), parameter :: PI = atan(1.0d0) * 4
  integer, parameter :: N = 16
  integer :: j, k

  integer(8) :: plan
  complex(8) :: in(0:N-1), out(0:N-1)

  ! 1. create a FFT plan first
  call dfftw_plan_dft_1d(plan, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE)

  do j = 0, N-1
    in(j) = sin(2 * PI * j / N)
  end do

  ! 2. execute FFT (as many times as you want)
  call dfftw_execute_dft(plan, in, out)

  do k = 0, N/2
    write(*, *) k, k, out(k)
  end do
  do k = N/2+1, N-1
    write(*, *) k, k - N, out(k)
  end do

  ! 3. destroy the plan finally
  call dfftw_destroy_plan(plan)

  stop
end program fftw_c2c
```

0. 最初にFFTについての**プラン**を作成する
	- 引数の `plan` はC言語のポインタ。64bit環境を考慮して `integer(8)` にしておく。
	- `N` の部分は入力配列の大きさとは限らない（後のc2rやDCT-Iなどが分かりやすい）。目的の変換をFFTで換算した場合の論理的な大きさを表す。
	- `in`, `out` にFFTの入出力配列を指定する。この時点ではデータを入れないこと。（内容が書き換えられる可能性がある）
	- 正変換（時間領域→周波数領域）では `FFTW_FORWARD`、逆変換（周波数領域→時間領域）では `FFTW_BACKWARD` 。
	- 最後の引数はFFTWの細かい設定をフラグの論理和（`ior()`）で指定する。フラグを二重に設定しない自信があれば単純に `+` でもいい。
0. そのプランを用いて、好きなだけ繰り返しFFTを実行する
	- 引数の配列はプラン作成時のものが安全。別の配列を用いるときは、必ずプラン作成時と同じ形状にする。
	- 入力配列の値が壊されることがあるので注意。（c2r, hc2rのデフォルトと、プランのフラグに `FFTW_DESTROY_INPUT` を指定して高効率と判断された場合）
0. 同プランのFFTを使わなくなったら破棄する

## 実数データのFFT
DFTは変換の前後でデータの自由度が変わらないため、入力値がすべて実数（虚部がゼロ）だと出力値の複素数は半分が冗長なデータになる。具体的には、**正波数と負波数の複素振幅は共役となる**（`out(modulo(-k, N)) == conjg(out(k))`）[^conjg]。この無駄な部分を削れば、計算時間とメモリ使用量を節約できる。

[^conjg]: 逆の波数が自分自身な `out(0)` などは虚部がゼロ

### r2c, c2r
複素数版から最も移行しやすい方法。時間領域のデータは実数配列になり、周波数領域のデータは**約半分（添字 0 ～ N/2 のみ）の複素数配列**になる。

プラン作成・実行時のサブルーチンが `dft` から `dft_r2c`, `dft_c2r` になり、入出力配列の引数の型がわかる。またこの名前自体が正変換・逆変換を表すので、`FFTW_FORWARD` のような種類の指定は引数に無い。

c2rはデフォルトで入力配列を書き換えるので、それが嫌なときはプラン作成のフラグに `FFTW_PRESERVE_INPUT` を加える。

### r2hc, hc2r
ここから先は**出力も実数配列**になる。プラン作成・実行時のサブルーチンは `dft` *ではなく* `r2r` で、種類の指定も復活する。

前項とは異なる方法で**実数配列に周波数成分を格納**する。並びに注意が必要だが、**cos成分とsin成分**に分解した（だから振幅が実数）と考えてもいい。

c2rと同様にhc2rはデフォルトで入力配列を書き換えるので、それが嫌なときはプラン作成のフラグに `FFTW_PRESERVE_INPUT` を加える。

### DCT, DST
入力する実数データに偶対称性/奇対称性がある場合、**出力はcos成分かsin成分だけ**になる。これらの場合はさらに最適化でき、半周期か1/4周期分だけの最小限なデータを用いた[離散コサイン変換](https://ja.wikipedia.org/wiki/離散コサイン変換)（DCT）・離散サイン変換（DST）がある。ただし、データ両端の偶奇性・半サンプル分のずれによって種類がたくさんある[^DCT_family]うえ、係数の流儀も多いので気を使わなければいけない。

FFTWを使う上では、**プラン作成時に指定する大きさが配列の大きさと異なる**点にも注意。基本的には配列長の2倍つまり 2N だが、DCT-Iなら 2(N-1)、DST-Iなら 2(N+1) となる。これらの値はFFTと関係していて、規格化係数であり、小さな素数の積のとき高速になる。



[^DCT_family]: 理論上は16通り、FFTWで使えるのは一般的な8通り

### DHT
実数データに対するFFT研究の歴史の中では、離散ハートレー変換（DHT）を使うほうが速い時期があったらしい。現在ではDFTに基づくほうが性能がいい[^FFTW_DHT]し、分解する波が少し異なるので、単純な性能目的で使う必要は無い。

規格化係数は配列長と同じ N 。正変換と逆変換が同じなので、種類は `FFTW_DHT` だけがある。

[^FFTW_DHT]: FFTWの内部では、r2hcで変換してから後処理でDHTに合わせている。（[マニュアル参照](http://www.fftw.org/fftw3_doc/The-Discrete-Hartley-Transform.html)）

## 雑記
本当は、畳み込み・微分・スペクトル計算への適用方法についても書きたかった（定数倍ずれる罠はそこにもある）。そこまで纏めきるのが難しかったため、今回はFFTWの使い方に留める。

<!--
.## 応用
定義では複素数のデータに対して計算できるが、実数のみを対象にすることで冗長な計算を省けることがある。よく分からなくなったら一旦複素数に立ち戻って冗長なところを確認すればいい。（複素数の計算式には複素共役が結構あるので注意）

.### [畳み込み](https://ja.wikipedia.org/wiki/畳み込み)
あるデータをずらしながら重ね合わせる演算。大雑把に言えば**巨大な掛け算の筆算**みたいなもの。重ね合わせの回数が少なければ時間領域でそのまま計算すればいいが、回数が多い（データの大きさと同じ）場合は畳み込み定理を利用して周波数領域での掛け算に持ち込む。

掛け算の筆算と違うのは、元の桁数をはみ出した部分が下位の桁に足されること（循環畳み込み）。その性質が好ましくない（線形畳み込みをしたい）場合は、データをゼロで拡張してから計算する。元のデータの2倍を超えていればいい。

FFTWを使う場合はどこかで規格化しないといけない。

```fortran
program sample
  implicit none
  include 'fftw3.f'

  integer, parameter :: N1 = 3, N2 = 3
  integer, parameter :: N = 8   ! N >= N1 + N2 - 1
  real(8)    :: rx(0:N-1), ry(0:N-1), rz(0:N-1)
  complex(8) :: cx(0:N/2), cy(0:N/2), cz(0:N/2)

  integer, parameter :: FLAGS = FFTW_ESTIMATE
  integer(8) :: plan_x, plan_y, plan_z

  call dfftw_plan_dft_r2c_1d(plan_x, N, rx, cx, FLAGS)
  call dfftw_plan_dft_r2c_1d(plan_y, N, ry, cy, FLAGS)
  call dfftw_plan_dft_c2r_1d(plan_z, N, cz, rz, FLAGS)

  ! (-1 + a + a^2) x (1 + a + a^2) = -1 + a^2 + 2a^3 + a^4
  rx(:) = 0.0d0;   rx(0:N1-1) = (/ -1.0d0, 1.0d0, 1.0d0 /)
  ry(:) = 0.0d0;   ry(0:N2-1) = (/  1.0d0, 1.0d0, 1.0d0 /)

  call dfftw_execute_dft_r2c(plan_x, rx, cx)
  call dfftw_execute_dft_r2c(plan_y, ry, cy)

  cz(:) = cx(:) * cy(:)

  call dfftw_execute_dft_c2r(plan_z, cz, rz)

  write(*, *) rz(0:N1+N2-2) / N

  call dfftw_destroy_plan(plan_x)
  call dfftw_destroy_plan(plan_y)
  call dfftw_destroy_plan(plan_z)

  stop
end program sample
```

.#### 周波数領域での畳み込み
逆に、周波数領域での畳み込みは時間領域での積になる。生成される高周波のエイリアシングを防ぐには同様にデータをゼロで拡張する必要があるが、もし高周波を切り落とすのならデータ数は元の 3/2 倍を超えていればいい（切り落とす周波数帯の中でならエイリアシングが起きてもいい）。

.### 微分
時間領域で微分するには近くのデータとの差分をとる必要があり、精度を上げるには複雑な計算が必要になる。しかし周波数領域に直すと振幅に周波数を掛けるだけで求まる。

```math
\frac{\partial}{\partial t} \left( A_\omega e^{i \omega t} \right)
= i \omega A_\omega e^{i \omega t}
```

実装ではデータの物理的な周期（`T = N * dt`）と正変換・逆変換の規格化に注意しないと倍率を間違えてしまう。

```fortran
program sample
  implicit none
  include 'fftw3.f'

  real(8), parameter :: PI = atan(1.0d0) * 4
  real(8), parameter :: T = 2 * PI * 5
  integer, parameter :: N = 16
  integer :: j, k

  integer(8) :: plan1, plan2
  real(8)    :: in(0:N-1)
  complex(8) :: out(0:N/2)

  call dfftw_plan_dft_r2c_1d(plan1, N, in, out, FFTW_ESTIMATE)
  call dfftw_plan_dft_c2r_1d(plan2, N, out, in, FFTW_ESTIMATE)

  do j = 0, N-1
    in(j) = sin(2 * PI * j / N)
  end do

  call dfftw_execute_dft_r2c(plan1, in, out)

  ! derivate
  do k = 0, N/2
    out(k) = out(k) * (0, 1) * k / N * (2 * PI / T)
  end do

  call dfftw_execute_dft_c2r(plan2, out, in)

  write(*, *) "# j, calc, exact"
  do j = 0, N-1
    write(*, *) j, in(j), cos(2 * PI * j / N) * (2 * PI / T)
  end do

  call dfftw_destroy_plan(plan1)
  call dfftw_destroy_plan(plan2)

  stop
end program sample
```

.### パワースペクトル密度
倍率が合っているかについては[パーセバルの定理](https://ja.wikipedia.org/wiki/パーセバルの定理)を使えば簡単に確認できる。

フーリエ変換周辺だけ実装しているので、データの前処理・後処理を適宜加えること。

.#### フーリエ変換法
データをフーリエ変換して周波数毎の振幅を直接手に入れる。FFTのおかげで高速にスペクトルが求まる。

.#### 自己相関法
[Wiener-Khinchinの定理](https://ja.wikipedia.org/wiki/ウィーナー＝ヒンチンの定理)を利用し、データの自己相関関数のフーリエ変換として求める。データに欠損が少量あっても計算できる。

実関数の自己相関関数は偶関数なので、DCT-Iを使えばメモリを節約できる。自己相関関数の計算に時間がかかるため、計算時間にはあまり効果がない。

```fortran
program sample
  implicit none
  include 'fftw3.f'

  real(8), parameter :: PI = atan(1.0d0) * 4
  real(8), parameter :: T = 2 * PI * 5
  integer, parameter :: N = 1024
  integer, parameter :: MAXLAG = N / 4   ! MAXLAG <= N / 2
  integer, parameter :: M = MAXLAG + 1
  integer :: j, k
  real(8) :: data(0:N-1)

  integer(8) :: plan
  real(8) :: corr(0:M-1), psd(0:M-1)

  call dfftw_plan_r2r_1d(plan, M, corr, psd, FFTW_REDEF00, FFTW_ESTIMATE)

  ! calc correlation function
  do k = 0, M-1
    corr(k) = sum(data(0:N-1-k) * data(k:N-1)) / (N - k)
  end do

  call dfftw_execute_r2r(plan, corr, psd)

  do k = 0, M-1
    write(*, *) k, psd(k)
  end do

  call dfftw_destroy_plan(plan)

  stop
end program sample
```
-->

## 参考

- [FFTW Manual](http://www.fftw.org/fftw3_doc/)
	- [Tutorial](http://www.fftw.org/fftw3_doc/Tutorial.html)
	- [What FFTW Really Computes](http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html)
	- [Calling FFTW from Legacy Fortran](http://www.fftw.org/fftw3_doc/Calling-FFTW-from-Legacy-Fortran.html)
- [FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法](http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/index.html)

<!-- - スペクトル法による数値計算入門 -->