scale_fast_fourier_transform.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> Module common / FFT
!!
!! @par Description
!!          A module to provide a fast Fourier transformation
!!
!! @author Yuta Kawai, Team SCALE
!!
!! @par 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. 
!<
#include "scaleFElib.h"
module scale_fast_fourier_transform
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_const, only:    &
    undef8 => const_undef8, &
    pi => const_pi
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
  type, public :: fastfouriertransform1d
    integer :: n
    integer :: m
    complex(RP), allocatable :: w_forward(:)
    complex(RP), allocatable :: w_backward(:)
 
    ! for N /= 2^n    
    logical :: use_bluestein
    complex(RP), allocatable :: w1_bluestein(:)
    complex(RP), allocatable :: w2_bluestein(:)
 
  contains
    procedure :: init => fastfouriertransform1d_init
    procedure :: final => fastfouriertransform1d_final
    procedure :: forward_cmplx => fastfouriertransform1d_forward_cmplx
    procedure :: forward_real => fastfouriertransform1d_forward_real
    generic :: forward => forward_cmplx, forward_real
    procedure :: backward => fastfouriertransform1d_backward
  end type fastfouriertransform1d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  
  !-----------------------------------------------------------------------------
contains
!OCL SERIAL
  subroutine fastfouriertransform1d_init( this, N )
    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
  end subroutine fastfouriertransform1d_init
 
!OCL SERIAL
  subroutine fastfouriertransform1d_forward_real( this, q, s )
    implicit none
    class(fastfouriertransform1d), intent(in) :: this
    real(RP), intent(in) :: q(this%N)
    complex(RP), intent(out) :: s(this%N)
 
    complex(RP) :: q_(this%N)
    integer :: i
    !-----------------------------------
 
    do i=1, this%N
        q_(i) = cmplx(q(i), 0.0_rp, kind=rp)
    end do
    call this%Forward_cmplx( q_, s )
    return
  end subroutine fastfouriertransform1d_forward_real
 
!OCL SERIAL
  subroutine fastfouriertransform1d_forward_cmplx( this, q, s )
    implicit none
    class(fastfouriertransform1d), intent(in) :: this
    complex(RP), intent(in) :: q(this%N)
    complex(RP), intent(out) :: s(this%N)
 
    integer :: i
    real(RP) :: scaling
    !-----------------------------------
 
    if ( this%use_bluestein ) then
      call fft1d_bluestein( s, &
        q, this%W_forward, this%W_backward, this%W1_bluestein, this%W2_bluestein, &
        this%N, this%M )
    else
      call fft1d_norecursive_core( s, &
        q, this%W_forward, this%N )
      ! call fft1d_recursive_core( s, &
      !   q, this%W_forward, this%N, this%N )
    end if
 
    scaling = 1.0_rp / real( this%N )
    s(:) = scaling * s(:)
    return
  end subroutine fastfouriertransform1d_forward_cmplx
 
!OCL SERIAL
  subroutine fastfouriertransform1d_backward( this, s, q )
    implicit none
    class(fastfouriertransform1d), intent(in) :: this
    complex(RP), intent(in) :: s(this%N)
    real(RP), intent(out) :: q(this%N)
 
    complex(RP) :: q_(this%N)
    integer :: i
    !-----------------------------------
 
    if ( this%use_bluestein ) then
      call fft1d_bluestein( q_, &
        s, this%W_forward, this%W_backward, this%W2_bluestein, this%W1_bluestein, &
        this%N, this%M )
    else
      call fft1d_norecursive_core( q_, &
        s, this%W_backward, this%N )
    end if  
 
    do i=1, this%N
        q(i) = real(q_(i), kind=rp)
    end do
    return
  end subroutine fastfouriertransform1d_backward
 
!OCL SERIAL
  subroutine fastfouriertransform1d_final( this )
    implicit none
    class(fastfouriertransform1d), intent(inout) :: this
    !-----------------------------------
 
    deallocate( this%W_forward, this%W_backward)
 
    if ( this%use_bluestein ) then
      deallocate( this%W1_bluestein, this%W2_bluestein )
    end if
    return
  end subroutine fastfouriertransform1d_final
 
!- private -------------
!OCL SERIAL
  subroutine fft1d_bluestein( x_out, x_in, W_forward, W_backward, W1, W2, N, M )
    implicit none
    integer, intent(in) :: N, M
    complex(RP), intent(out) :: x_out(N)
    complex(RP), intent(in) :: x_in(N)
    complex(RP), intent(in) :: W_forward(M/2)  !< twiddle factor for forward transformation
    complex(RP), intent(in) :: W_backward(M/2) !< twiddle factor for backward transformation
    complex(RP), intent(in) :: W1(N)
    complex(RP), intent(in) :: W2(N)
 
    complex(RP) :: a(M), b(M)
    complex(RP) :: AA(M), BB(M), CC(M)
    integer :: i
    real(RP) :: scaling
    !---------------------------
 
    !- Chirp pre-processing
    a(:) = 0.0_rp
    b(:) = 0.0_rp
    do i=0, n-1
      a(i+1) = x_in(i+1) * w1(i+1)
      b(i+1) = w2(i+1)
    end do
    ! Fill symmetric part of b
    do i=1, n-1
      b(m-i+1) = w2(i+1)
    end do
 
    !- Convolution C = A * B 
 
    call fft1d_norecursive_core( aa, &
      a, w_forward, m )
    call fft1d_norecursive_core( bb, &
      b, w_forward, m )
 
    call fft1d_norecursive_core( cc, &
      aa(:) * bb(:), w_backward, m )
    
    !- Chirp post-processing
    scaling = 1.0_rp / real( m )    
    do i=1, n
      x_out(i) = scaling * cc(i) * w1(i)
    end do
    return
  end subroutine fft1d_bluestein
 
!OCL SERIAL
  subroutine fft1d_norecursive_core( x_out, x_in, W, N )
    implicit none
    integer, intent(in) :: N
    complex(RP), intent(out) :: x_out(N)
    complex(RP), intent(in) :: x_in(N)
    complex(RP), intent(in) :: W(N/2) !< twiddle factor
 
    integer :: rev
    integer :: nbits
    complex(RP) :: tmp(N)
 
    integer :: i, j
    integer :: k
    integer :: m
    integer :: ind
 
    integer :: stage
    integer :: step
 
    integer :: s
    complex(RP) :: t, u
    !---------------------------
 
    !- bit-reversal permutation
    tmp(:) = x_in(:)
 
    nbits = int( log(real(n,kind=rp)) / log(2.0_rp) + 0.5_rp )
    do i=0, n-1
      rev = 0
      do j=0, nbits-1
        if ( i/2**j - 2*(i/2**(j+1)) == 1 ) rev = rev + 2**(nbits-j-1)
      end do
      x_out(rev+1) = tmp(i+1)
    end do
 
    !- Iterative butterfly using W
    m = 1
    do stage=1, nbits
      step = 2*m
      s = n / step
      do k=0, n-1, step
        do j=0, m-1
          ind = j*s
          t = w(ind+1) * x_out(k+j+m+1)
          u = x_out(k+j+1)
          x_out(k+j+1  ) = u + t
          x_out(k+j+1+m) = u - t
        end do
      end do
      m = step
    end do
    return
  end subroutine fft1d_norecursive_core
 
!OCL SERIAL
  recursive subroutine fft1d_recursive_core( x_out, x_in, W, N, N0 )
    implicit none
    integer, intent(in) :: N
    integer, intent(in) :: N0
    complex(RP), intent(out) :: x_out(N)
    complex(RP), intent(in) :: x_in(N)
    complex(RP), intent(in) :: W(N0/2) !< twiddle factor
 
    complex(RP) :: even(N/2), odd(N/2)
    complex(RP) :: E(N/2), O(N/2)
 
    integer :: k
    integer :: s
    complex(RP) :: twiddle
    !---------------------------
 
    if ( n==1 ) then
      x_out(1) = x_in(1)
      return
    end if
 
    even(:) = x_in(1:n:2)
    odd(:) = x_in(2:n:2)
 
    call fft1d_recursive_core( e, &
      even, w, n/2, n0 )
    call fft1d_recursive_core( o, &
       odd, w, n/2, n0 )
 
    s = n0 / n
    do k=0, n/2-1
      x_out(k+1    ) = e(k+1) + w(s*k+1) * o(k+1)
      x_out(k+1+n/2) = e(k+1) - w(s*k+1) * o(k+1)
    end do
    return
  end subroutine fft1d_recursive_core
 
!OCL SERIAL
  logical function is_power_of_two(N)
    implicit none
    integer, intent(in) :: N
    !--------------
    is_power_of_two = (n > 0) .and. (iand(n, n-1) == 0)
    return
  end function is_power_of_two  
end module scale_fast_fourier_transform