FFTW3 でフーリエ変換

インストール

本家HP (リンク)に移動して http: fftw-3.3.10.tar.gz をダウンロード、適当な場所に移動しておく。

tar -zxvf http: fftw-3.3.10.tar.gz 
cd fftw-3.3.10
./configure --prefix=/usr/local/
make
sudo make install

他にも並列化オプションなどが必要であれば設定しておく。

ソースコード・コンパイル

10Hz, 20Hz, 30Hzの sin 信号を足し合わせた以下の信号をFFTする:
$$f(t) = \sin(2\pi 10t) + 0.7\sin(2\pi 20t) + 0.3\sin(2\pi 30t)$$ もちろん、FFTスペクトルにはf=10, 20, 30 に y= 1, 0.7, 0.3 のピークが立つ.

program test
implicit none

include "fftw3.f"

integer, parameter :: N = 100
real(8), parameter :: dt=1.d0/dble(N-1), fs=1.d0/dt
real(8)    :: t(N), freq(N), f(N)
complex(8) :: cf(N)     ! f(t) を複素数配列に変換したもの
complex(8) :: f1(N)     ! FFT[ f(t) ]
complex(8) :: f2(N)     ! INVFFT[ FFT[ f(t) ] ]
integer(8) :: plan
real(8), parameter :: K=2.d0*4.d0*atan(1.d0)
integer :: i

! set f(t)
do i = 1, N
    t(i)    = real(i-1) * dt
    freq(i) = dble(i-1) * fs / dble(N-1)
    f(i)    = sin(K*10.d0*t(i)) + 0.7d0*sin(K*20.d0*t(i)) + 0.3d0*sin(K*30.d0*t(i))
    cf(i)   = cmplx(f(i), 0.0d0) ! real -->> complex
end do

! FFT[ f(t) ]
call dfftw_plan_dft_1d(plan, N, cf, f1, FFTW_FORWARD, FFTW_ESTIMATE)
call dfftw_execute_dft(plan, cf, f1)
call dfftw_destroy_plan(plan)
call dfftw_cleanup()

! inversion FFT
call dfftw_plan_dft_1d(plan, N, f1, f2, FFTW_BACKWARD, FFTW_ESTIMATE)
call dfftw_execute_dft(plan, f1, f2)
call dfftw_destroy_plan(plan)
call dfftw_cleanup()

! amplitude setting
f1 = f1/dble(N/2)
f2 = f2/dble(N)

do i = 1, N
    write(55,'(3ES16.2)') t(i), f(i), real(f2(i))
    write(56,'(4ES16.2)') freq(i), real(f1(i)), aimag(f1(i)), abs(f1(i))
end do

end program test

gfortran test.f90 -I/usr/local/include/ -L/usr/local/library -lfftw3
./a.out

結果

問題なく計算できた。

エイリアシングも見えている。

その他

program test
implicie none
include "fftw3.f"

のように fftw3.f (変数定義のみされたファイル) を include する方式は使いにくい。
できれば use fftw3 のようにしたいところ。
この場合、fftw3.f の先頭行に module fftw3、末行に end module fftw3 と書いて、別途コンパイルしておいてもよいかも。