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の種類に応じて、サブルーチン名や引数、配列の型と大きさがどうなっていればいいかを纏める。
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のところで指定する
$ 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
$ ls ${HOME}/local/include #=> fftw3.f
$ ls ${HOME}/local/lib #=> libfftw3.a
コンパイルとリンク
- プログラムにはFortran用のincludeファイル
fftw3.fを読み込ませる- 拡張子からわかる通りFORTRAN77形式だが、Fortran90以降でも問題なく使える
- コンパイル時にincludeファイルの場所を教える
- リンク時にライブラリ
libfftw3.aの場所を教え、使用を明示する- コンパイルとリンクを同時にする場合は両方の指定が必要
- 環境変数の設定によっては
-I,-Lは不要
program sample
implicit none
include 'fftw3.f'
write(*, *) FFTW_ESTIMATE
stop
end program sample
$ 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
最も基本的な使い方。効率を気にしないならこれだけでも行ける。
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))）1。この無駄な部分を削れば、計算時間とメモリ使用量を節約できる。
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）がある。ただし、データ両端の偶奇性・半サンプル分のずれによって種類がたくさんある2うえ、係数の流儀も多いので気を使わなければいけない。
FFTWを使う上では、プラン作成時に指定する大きさが配列の大きさと異なる点にも注意。基本的には配列長の2倍つまり 2N だが、DCT-Iなら 2(N-1)、DST-Iなら 2(N+1) となる。これらの値はFFTと関係していて、規格化係数であり、小さな素数の積のとき高速になる。
DHT
実数データに対するFFT研究の歴史の中では、離散ハートレー変換（DHT）を使うほうが速い時期があったらしい。現在ではDFTに基づくほうが性能がいい3し、分解する波が少し異なるので、単純な性能目的で使う必要は無い。
規格化係数は配列長と同じ N 。正変換と逆変換が同じなので、種類は FFTW_DHT だけがある。
雑記
本当は、畳み込み・微分・スペクトル計算への適用方法についても書きたかった（定数倍ずれる罠はそこにもある）。そこまで纏めきるのが難しかったため、今回はFFTWの使い方に留める。