FFTWをfortranで使う。
本家の説明:
FFTW 3.3.8: Fortran Examples
参考にしたサイト:
fftwの使い方:腰も砕けよ 膝も折れよ:So-net blog
FortranのFFTWの使い方 - Qiita
Hiroyuki R. Takahashi
Fumio KUSUHARA -Fortran90で構造解析-
program fftw3_test implicit none integer, parameter :: WP = selected_real_kind( 8 ) real( WP ), parameter :: PI = 3.141592653589793 integer :: i1 ! setting parameters and variables integer, parameter :: N = 2**10 real( WP ) , parameter :: L = 100. real( WP ) :: dx real( WP ) :: x( N ), k( N ) complex( WP ) :: diff ! data complex( WP ) :: input0( N ), input( N ), output( N ) ! plan integer( WP ) :: plan include 'fftw3.f' dx = L / N ! test: Gaussian function do i1 = 1, N x( i1 ) = - L/2 + ( i1 - 1 ) * dx input( i1 ) = exp( - x( i1 )**2 ) enddo input0( : ) = input( : ) call dfftw_plan_dft_1d( plan, N, input0, output, & & FFTW_FORWARD, FFTW_ESTIMATE ) call dfftw_execute_dft( plan, input0, output ) call dfftw_destroy_plan( plan ) do i1 = 1, N/2 k( i1 ) = i1 - 1 k( N - ( i1 - 1 ) ) = - i1 enddo k( : ) = ( 2 * PI / L ) * k( : ) open( 10, file = 'fftw3_forward.out', status = 'unknown' ) do i1 = 1, N write( 10, '(4e15.5)' ) k( i1 ), & & real( output( i1 ) ), aimag( output( i1 ) ), abs( output( i1 ) ) enddo close( 10 ) call dfftw_plan_dft_1d( plan, N, output, input, & & FFTW_BACKWARD, FFTW_ESTIMATE ) call dfftw_execute_dft( plan, output, input ) call dfftw_destroy_plan( plan ) input( : ) = input( : ) / N open( 10, file = 'fftw3_backward.out', status = 'unknown' ) do i1 = 1, N diff = input( i1 ) - input0( i1 ) write( 10, '(4e15.5)' ) x( i1 ), & & real( diff ), aimag( diff ), abs( diff ) enddo close( 10 ) end program fftw3_test
フーリエ変換後の配列は、0, 1, 2, ..., N/2, -(N/2 - 1), -(N/2 - 2), ... -1 となっていることに注意。
この配列のまま逆フーリエ変換にかけると元に戻る。
compile:
gfortran -o fftw3_test.out fftw3_test.f90 -I/usr/local/include -L/usr/local/lib -lfftw3