Below is my wrapper for the ooura real discrete Fourier transform. The data to be transformed, including real and complex pointers to the same data, are components of rfft_type from which rfft_ooura_type extends. According to gfortran I am not allowed to have allocatable targets as type components so I had to allocate the cview pointer (against FortranFan’s advice…). Note that size(cview) equals size(rview)/2+1. The reason for the extra element is that the imaginary part of the zero and Nyquist frequencies are always zero. The rdft routine of Ooura puts the real part of the Nyquist component in the position of the imaginary part of the zero frequency. My wrapper moves the Nyquist component to the end where it belongs, thereby behaving like the rfft of Python-numpy.
I could have made rview and cview private components with getter and setter methods, but that will probably lead to copying in and out(?) which is what I want to avoid in the first place…
subroutine test_rfft_ooura()
type(rfft_ooura_type) :: fobj
integer nt, i, ierr
nt = 4
call fobj%initialize(nt, ierr)
! Create time series
fobj%rview = [ (real(i, wp) , i=1,nt)]
print*,'Time series with size ',size(fobj%rview)
print*,'rview: ',fobj%rview
! Transform to frequency domain
call fobj%rfft()
print*,'Fourier transform with size ',size(fobj%cview)
print*,'cview: ',fobj%cview
! Transform back to time domain
call fobj%irfft()
print*,'Time series (recreated) '
print*,'rview: ',fobj%rview
call fobj%destroy()
end subroutine test_rfft_ooura
Output
Time series with size 4
rview: 1.0000000000000000 2.0000000000000000 3.0000000000000000 4.0000000000000000
Fourier transform with size 3
cview: (10.000000000000000,0.0000000000000000) (-2.0000000000000000,2.0000000000000000) (-2.0000000000000000,0.0000000000000000)
Time series (recreated)
rview: 1.0000000000000000 2.0000000000000000 3.0000000000000000 4.0000000000000000
module rfft_ooura
use fftsg_modified_mod, only : rdft
use rfft, only : rfft_type
implicit none
private
public :: rfft_ooura_type , initialize, destroy, is_initialized, rfft, irfft
integer, parameter :: wp=kind(0.d0) !Move this to subroutines
type, extends(rfft_type) :: rfft_ooura_type
private
logical :: initialized = .false.
integer :: nr, nc
integer, allocatable :: ip(:)
real(wp), allocatable :: cstab(:)
! From parent: complex(wp), public, pointer, contiguous :: cview(:)
! From parent: real(wp), public, pointer :: rview(:)
contains
procedure :: initialize, destroy, is_initialized, rfft, irfft
end type rfft_ooura_type
contains
subroutine initialize(this, fftsize, ierr)
use ISO_C_binding, only: c_f_pointer, c_loc
use fft_utils_mod, only: next_power_of_2
class(rfft_ooura_type), intent(inout) :: this
integer, intent(in) :: fftsize
integer, intent(out) :: ierr
if (fftsize /= next_power_of_2(fftsize)) then
ierr = -1
return
else
ierr = 0
end if
! Set status
this%initialized = .true.
this%nr = fftsize
this%nc = fftsize/2 + 1
! Allocate work arrays
allocate(this%ip(2+2**(int(log(fftsize/2+0.5)/log(2.0))/2)))
allocate(this%cstab(fftsize/2))
this%ip(1) = 0 !Initialization flag
! Allocate data (to be transformed) with both a real and a complex "view"
! (rview and cview is pointers to the same data)
allocate(this%cview(this%nc))
call c_f_pointer(c_loc(this%cview), this%rview, shape=[fftsize])
this%cview = 0
end subroutine
!*****************************************************************************
subroutine destroy(this)
class(rfft_ooura_type), intent(inout) :: this
this%initialized = .false.
deallocate(this%ip, this%cstab, this%cview)
nullify(this%cview, this%rview)
end subroutine
!*****************************************************************************
function is_initialized(this)
class(rfft_ooura_type), intent(in) :: this
logical :: is_initialized
is_initialized = this%initialized
end function
!*****************************************************************************
subroutine rfft(this)
class(rfft_ooura_type), intent(inout) :: this
call rdft(this%nr, 1, this%rview, this%ip, this%cstab)
! Move the Nyquist component where it should be:
this%cview(this%nc)%re = this%cview(1)%im
this%cview(this%nc)%im = 0
this%cview(1)%im = 0
! Flip sign of imaginary part
this%cview(2:(this%nc-1))%im = -this%cview(2:(this%nc-1))%im
end subroutine
!*****************************************************************************
subroutine irfft(this)
class(rfft_ooura_type), intent(inout) :: this
! Flip sign of imaginary part
this%cview(2:(this%nc-1))%im = -this%cview(2:(this%nc-1))%im
! Move the Nyquist component where rdft expects it to be:
this%cview(1)%im = this%cview(this%nc)%re
call rdft(this%nr, -1, this%rview, this%ip, this%cstab)
this%rview = this%rview*(2.0_wp/this%nr) !scaling
this%cview(this%nc) = 0 ! Set extra element to zero by definition
end subroutine
!*****************************************************************************
end module rfft_ooura
Here is the abstract type:
module rfft
implicit none
private
public :: rfft_type
integer, parameter :: wp=kind(0.d0) !Move this to subroutines
type, abstract :: rfft_type
complex(wp), public, pointer, contiguous :: cview(:)
real(wp), public, pointer :: rview(:)
contains
private
! Initialization may vary between subtypes and is therefore not here
procedure(destroy_abstract), deferred :: destroy
procedure(is_initialized_abstract), deferred :: is_initialized
procedure(rfft_abstract), deferred :: rfft
procedure(irfft_abstract), deferred :: irfft
end type rfft_type
interface
subroutine destroy_abstract(this)
import :: rfft_type
class(rfft_type), intent(inout) :: this
end subroutine
end interface
interface
function is_initialized_abstract(this) result(is_initialized)
import :: rfft_type
class(rfft_type), intent(in) :: this
logical :: is_initialized
end function
end interface
interface
subroutine rfft_abstract(this)
import :: rfft_type
class(rfft_type), intent(inout) :: this
end subroutine
end interface
interface
subroutine irfft_abstract(this)
import :: rfft_type
class(rfft_type), intent(inout) :: this
end subroutine
end interface
end module rfft