Fortran 95 source code to compute discrete Fourier transformation.

discrete_fourier_transformation.f95

!****************************************************
! 離散フーリエ変換

! * 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