FFTWは、高速フーリエ変換（FFT）に基づいて離散フーリエ変換（DFT）を高速に計算できるライブラリのひとつ。任意の長さのデータが扱え、多次元FFTや実数に最適化した変換も用意されている。

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

実行環境など

記事では以下の環境を用いている。

FFTW 3.3.8
言語: Fortran90 ← "Calling FFTW from Legacy Fortran"
コンパイラ: gcc (gfortran含む)
OS: Ubuntu 16.04

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

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

TL;DR: プラン作成一覧

FFTの種類に応じて、サブルーチン名や引数、配列の型と大きさがどうなっていればいいかを纏める。

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

FFTWの使用準備

既にコンパイルまでできるなら飛ばしていい。

インストール

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

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

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

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 は不要

sample.f90

program sample
  implicit none
  include 'fftw3.f'

  write(*, *) FFTW_ESTIMATE

  stop
end program sample

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

最も基本的な使い方。効率を気にしないならこれだけでも行ける。

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

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

実数データのFFT

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

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周期分だけの最小限なデータを用いた離散コサイン変換（DCT）・離散サイン変換（DST）がある。ただし、データ両端の偶奇性・半サンプル分のずれによって種類がたくさんある²うえ、係数の流儀も多いので気を使わなければいけない。

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

DHT

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

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

雑記

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

参考

逆の波数が自分自身な out(0) などは虚部がゼロ ↩
理論上は16通り、FFTWで使えるのは一般的な8通り ↩
FFTWの内部では、r2hcで変換してから後処理でDHTに合わせている。（マニュアル参照） ↩