scale_fast_fourier_transform Module Reference

Module common / FFT. More...

Data Types
type	fastfouriertransform1d

Functions/Subroutines
subroutine	fastfouriertransform1d_init (this, n)

Detailed Description

Module common / FFT.

Description

A module to provide a fast Fourier transformation

Author

Yuta Kawai, Team SCALE

Reference

• Cooleyand & Tukey, 1965: An algorithm for the machine calculation of complex Fourier series Mathematics of Computation, 19, 297-301
• Bluestein 1970: A linear filtering approach to the computation of discrete Fourier transform IEEE Transactions on Audio and Electroacoustics, 18(4), 451-455.

Function/Subroutine Documentation

fastfouriertransform1d_init()

subroutine scale_fast_fourier_transform::fastfouriertransform1d_init	(	class(fastfouriertransform1d), intent(inout)	this,
		integer, intent(in)	n )

Definition at line 74 of file scale_fast_fourier_transform.F90.

    implicit none
    class(FastFourierTransform1D), intent(inout) :: this
    integer, intent(in) :: N
 
    complex(RP) :: JJ
    integer :: k
 
    real(RP) :: angle0
    real(RP) :: angle
 
    integer :: M
    integer :: MM
 
    complex(RP) :: w1, w2
    !-----------------------------------
 
    jj = cmplx( 0.0_rp, 1.0_rp )
 
    this%N = n
    if ( is_power_of_two(n) ) then
      this%use_bluestein = .false.
      mm = n / 2
      log_info("FastFourierTransform1D_Init",*) "use_bluestein=F"
      log_info("FastFourierTransform1D_Init",*) "MM=", mm
    else
      this%use_bluestein = .true.
      m = 1
      do while ( m < 2*n-1 )
        m = m * 2
      end do
      mm = m / 2
      log_info("FastFourierTransform1D_Init",*) "use_bluestein=T"
      log_info("FastFourierTransform1D_Init",*) "M=", m, "MM=", mm
    end if
    this%M = m
 
    allocate( this%W_forward(mm), this%W_backward(mm) )
 
    angle0 = 2.0_rp * pi / real(2*mm, kind=rp)
    !$omp parallel do private(angle)
    do k=0, mm-1
      angle = angle0 * real(k, kind=rp)
      this%W_forward (k+1) = exp(- jj * angle )
      this%W_backward(k+1) = exp(  jj * angle )
    end do
 
    if ( this%use_bluestein ) then
      allocate( this%W1_bluestein(n), this%W2_bluestein(n) )
      !$omp parallel do private(angle)
      do k=0, n-1
        angle = pi * real(k*k,kind=rp) / real(n, kind=rp)
        this%W1_bluestein(k+1) = exp(- jj * angle )
        this%W2_bluestein(k+1) = exp(  jj * angle )
      end do
    end if
 
    return