fortran90で高速フーリエ変換

Program FFT

! FFT(高速フーリエ変換)
!くぴん７４
! 2012.9.27 移植
! 2014.6.2 移植
! 2014.6.17 完成
! 2014.7.5 修正

implicit none
double precision cr(0:100, 0:10000), ci(0:100, 0:10000)
double precision co(0:10000), si(0:10000)
double precision pi
integer n(0:10000)
integer l !ポイント数(2の乗数)
integer m
double precision:: fs=1000. !サンプリング周波数[Hz]
double precision log2

pi = 4.*atan(1.) !円周率定義
Call data_input(cr, ci, l, m) !データ入力
Call bit_reverse1(co, si, pi, l, n) !nのビットリバースとωの作成
Call butterfly(cr, ci, co, si, l, m, n) !バタフライ演算
Call bit_reverse2(cr, ci, l, n) !スペクトルのビットリバース
Call spectrum(cr, ci, pi, l, fs) !スペクトルの出力

End Program FFT

!--------------------- データ入力 --------------------------
subroutine data_input(cr, ci, l, m)

implicit none
double precision cr(0:100, 0:10000), ci(0:100, 0:10000)
integer i
integer l !ポイント数(2の乗数)
integer l0
integer m
integer io
double precision:: fs=1000. !サンプリング周波数[Hz]
double precision log2
character(1000) file_in
! file_in="/home/pyontaku14/fortran_practice/10_50_100Hz.txt"
! 引数で入力ファイルを取得
call getarg(1, file_in)
open(3, file=file_in, status='old', err=99)
do i=0,10000
read(3,*, iostat=io, end=100) cr(0, i)
end do
99 write(*,*) "no file"
100 write(*,*) "end reading", i
close(3)
l = i
m = l
i = 0
do while (2 ** i < l)
i = i + 1
end do
l0 = l
l = 2 ** i
m = l
do i = l0 + 1, l
cr(0, i) = 0.0d0
end do

End subroutine
!----------------- nのビットリバースとωの作成 ---------------
subroutine bit_reverse1(co, si, pi, l, n)

implicit none
double precision co(0:10000), si(0:10000)
double precision pi
integer n(0:10000)
integer l !ポイント数(2の乗数)
integer j, k
double precision log2

do k = 0, l - 1
co(k) = cos(2.0d0 * dble(k) * pi / dble(l))
si(k) = sin(2.0d0 * dble(k) * pi / dble(l))
do j = 0, nint(log2(dble(l)) - 1.0d0)
n(k) = n(k) + mod((k / 2 ** j), 2) * nint(2.0d0 ** (log2(dble(l)) - 1.0d0 - dble(j)))
end do
end do

End subroutine
!---------------------- バタフライ演算 ----------------------
subroutine butterfly(cr, ci, co, si, l, m, n)

implicit none
double precision cr(0:100, 0:10000), ci(0:100, 0:10000)
double precision co(0:10000), si(0:10000)
integer n(0:10000)
integer l !ポイント数(2の乗数)
integer m
integer j, jj, k, kk
double precision log2

do j = 0, nint(log2(dble(l)) - 1.0d0)
kk = 2 ** j - 1
m = m / 2
do jj = 0, kk
do k = 0, m - 1
cr(j + 1, k + 2 * jj * m) = cr(j, k + 2 * jj * m) &
+ cr(j, k + 2 * jj * m + m) * co(n(2 * jj)) + ci(j, k + 2 * jj * m + m) * si(n(2 * jj))
ci(j + 1, k + 2 * jj * m) = ci(j, k + 2 * jj * m) - cr(j, k + 2 * jj * m + m) * si(n(2 * jj)) &
+ ci(j, k + 2 * jj * m + m) * co(n(2 * jj))
cr(j + 1, k + 2 * jj * m + m) = cr(j, k + 2 * jj * m) &
+ cr(j, k + 2 * jj * m + m) * co(n(2 * jj + 1)) + ci(j, k + 2 * jj * m + m) * si(n(2 * jj + 1))
ci(j + 1, k + 2 * jj * m + m) = ci(j, k + 2 * jj * m) &
- cr(j, k + 2 * jj * m + m) * si(n(2 * jj + 1)) + ci(j, k + 2 * jj * m + m) * co(n(2 * jj + 1))
end do
end do
end do

End subroutine
!----------------- スペクトルのビットリバース ----------------
subroutine bit_reverse2(cr, ci, l, n)

implicit none
double precision cr(0:100, 0:10000), ci(0:100, 0:10000)
integer n(0:10000)
integer l !ポイント数(2の乗数)
integer k
double precision log2

do k = 0, l - 1
cr(nint(log2(dble(l)) + 1.0d0), n(k)) = cr(nint(log2(dble(l))), k)
ci(nint(log2(dble(l)) + 1.0d0), n(k)) = ci(nint(log2(dble(l))), k)
end do

End subroutine
!--------------------- スペクトルの出力 ----------------------
subroutine spectrum(cr, ci, pi, l, fs)

implicit none
double precision cr(0:100, 0:10000), ci(0:100, 0:10000)
double precision pi
integer l !ポイント数(2の乗数)
integer k
double precision fs !サンプリング周波数[Hz]
double precision log2
character(1000) file_out

! 出力ファイル
file_out="fft_out.txt"
open(4, file=file_out, status='replace')

do k = 0, l - 1
! 周波数、実部、虚部、振幅、位相
write(4,"(f15.6, f15.6,f15.6,f15.6,f15.6)") &
dble(k) * fs / dble(l) , &
(2.0d0 / dble(l)) * cr(nint(log2(dble(l)) + 1.0d0), k), &
(2.0d0 / dble(l)) * ci(nint(log2(dble(l)) + 1.0d0), k), &
(2.0d0 / dble(l)) * sqrt(cr(nint(log2(dble(l)) + 1.0d0), k) ** 2 + ci(nint(log2(dble(l)) + 1.0d0), k) ** 2), &
(180.0d0 / pi) * atan2(ci(nint(log2(dble(l)) + 1.0d0), k), cr(nint(log2(dble(l)) + 1.0d0), k))
end do
close(4)

End subroutine

function log2(x)
implicit none
double precision x, log2
log2=log(x)/log(2.0d0)
end function