Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to compute discrete Fourier transformation.
!****************************************************
! 離散フーリエ変換
! * f(t) = 2 * sin(4 * t) + 3 * cos(2 * t)
! ( 0 <= t < 2 * pi )
!
! date name version
! 2018.12.17 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
real(DP), parameter :: PI = 4.0_DP * atan(1.0_DP) ! 円周率
integer(SP), parameter :: NUM = 100 ! 分割数
end module const
module fourier
use const
implicit none
private
public :: gen_src, exec_dft, exec_idft
contains
! 元データ生成
!
! :param(in) integer(4) n: データ個数(行数)
! :param(inout) real(8) dat(7, n): データ配列
subroutine gen_src(n, dat)
implicit none
integer(SP), intent(in) :: n
real(DP), intent(inout) :: dat(7, 0:n-1)
integer(SP) :: i
do i = 0, n - 1
dat(1, i) = (2 * PI / n) * i
dat(2, i) = 2 * sin(4 * (2 * PI / n) * i) &
& + 3 * cos(2 * (2 * PI / n) * i)
dat(3, i) = 0.0_DP
end do
end subroutine gen_src
! 離散フーリエ変換
!
! :param(in) integer(4) n: データ個数(行数)
! :param(inout) real(8) dat(7, n): データ配列
subroutine exec_dft(n, dat)
implicit none
integer(SP), intent(in) :: n
real(DP), intent(inout) :: dat(7, 0:n-1)
integer(SP) :: i, j
real(DP) :: dft_re, dft_im ! 変換後(実部、虚部)
do i = 0, n - 1
dft_re = 0.0_DP
dft_im = 0.0_DP
do j = 0, n - 1
dft_re = dft_re &
& + dat(2, j) * ( cos((2 * PI / n) * i * j)) &
& + dat(3, j) * ( sin((2 * PI / n) * i * j))
dft_im = dft_im &
& + dat(2, j) * (-sin((2 * PI / n) * i * j)) &
& + dat(3, j) * ( cos((2 * PI / n) * i * j))
end do
dat(4, i) = dft_re
dat(5, i) = dft_im
end do
end subroutine exec_dft
! 逆離散フーリエ変換
!
! :param(in) integer(4) n: データ個数(行数)
! :param(inout) real(8) dat(7, n): データ配列
subroutine exec_idft(n, dat)
implicit none
integer(SP), intent(in) :: n
real(DP), intent(inout) :: dat(7, 0:n-1)
integer(SP) :: i, j
real(DP) :: idft_re, idft_im ! 逆変換後(実部、虚部)
do j = 0, n - 1
idft_re = 0.0_DP
idft_im = 0.0_DP
do i = 0, n - 1
idft_re = idft_re &
& + dat(4, i) * (cos((2 * PI / n) * i * j)) &
& - dat(5, i) * (sin((2 * PI / n) * i * j))
idft_im = idft_im &
& + dat(4, i) * (sin((2 * PI / n) * i * j)) &
& + dat(5, i) * (cos((2 * PI / n) * i * j))
end do
dat(6, j) = idft_re / n
dat(7, j) = idft_im / n
end do
end subroutine exec_idft
end module fourier
program discrete_fourier_transformation
use const, only : SP, DP, NUM
use fourier
implicit none
real(DP) :: dat(7, 0:NUM-1) ! 計算データ用配列
integer(SP) :: i
dat(:, :) = 0.0_DP ! 配列初期化
call gen_src(NUM, dat) ! 元データ作成
call exec_dft(NUM, dat) ! 離散フーリエ変換
call exec_idft(NUM, dat) ! 逆離散フーリエ変換
! 結果出力
! * 左から:f, 元(実), (虚), IFT(実), (虚), DIFT(実), (虚)
do i = 0, NUM - 1
print '(7(X, F11.6))', dat(:, i)
end do
end program discrete_fourier_transformation
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.