module fft_work use dcl_common type work integer :: n real, dimension(:), pointer :: array end type work end module !------------------------------------------------- ! FFT 周期実数値データのフーリエ変換 (1) !------------------------------------------------- module fftreal use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitRealFFT(n, index) !初期化をおこなう integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitRealFFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 2*n+15 allocate(wk(idx)%array(len)) call rffti(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocRealFFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- function DclRealFFT_F(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclRealFFT_F if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclRealFFT_F', & & 'Wrong working area.') DclRealFFT_F = r call rfftf(n, DclRealFFT_F, wk(idx)%array) end function !--------------------------------------------------------- function DclRealFFT_B(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclRealFFT_B if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclRealFFT_B', & & 'Wrong working area.') DclRealFFT_B = r call rfftb(n, DclRealFFT_B, wk(idx)%array) end function end module !------------------------------------------------- ! FFT: rffti, rfftf, rfftbの簡易型サブルーチン (2) !------------------------------------------------- module ffteasy use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitEasyFFT(n, index) !初期化をおこなう integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitEasyFFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 3*n+15 allocate(wk(idx)%array(len)) call ezffti(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocEasyFFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclEasyFFT_F(r,a0,a,b,index) real, intent(in), dimension(:) :: r real, intent(out) :: a0 real, intent(out), dimension(:) :: a real, intent(out), dimension(:) :: b integer, intent(in), optional :: index if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclEasyFFT_F', & & 'Wrong working area.') call ezfftf(n,r,a0,a,b,wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclEasyFFT_B(r,a0,a,b,index) real, intent(out), dimension(:) :: r real, intent(in) :: a0 real, intent(in), dimension(:) :: a real, intent(in), dimension(:) :: b integer, intent(in), optional :: index if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclEasyFFT_B', & & 'Wrong working area.') call ezfftb(n,r,a0,a,b,wk(idx)%array) end subroutine end module !------------------------------------------------- ! 奇の周期データのsine変換をおこなう. (3) !------------------------------------------------- module fftsin use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitSinFFT(n,index) integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitSinFFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 2.5*n+15 allocate(wk(idx)%array(len)) call sinti(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocSinFFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- function DclSinFFT(r,index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclSinFFT if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinFFT', & & 'Wrong working area.') DclSinFFT = r call sint(n, DclSinFFT, wk(idx)%array) end function end module !------------------------------------------------- ! FFT: 偶の周期データのcosine変換をおこなう (4) !------------------------------------------------- module fftcos use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitCosFFT(n,index) integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitCosFFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 3*n+15 allocate(wk(idx)%array(len)) call costi(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocCosFFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- function DclCosFFT(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclCosFFT if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosFFT', & & 'Wrong working area.') DclCosFFT = r call cost(n,DclCosFFT,wk(idx)%array) end function end module !------------------------------------------------- ! FFT: 奇数波数成分のみのsin変換をおこなう.(5) !------------------------------------------------- module fftqsin use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitSinQFT(n,index) integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitSinQFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 3*n+15 allocate(wk(idx)%array(len)) call sinqi(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocSinQFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- function DclSinQFT_F(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclSinQFT_F if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinQFT_F', & & 'Wrong working area.') DclSinQFT_F = r call sinqf(n, DclSinQFT_F, wk(idx)%array) end function !--------------------------------------------------------- function DclSinQFT_B(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclSinQFT_B if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinQFT_B', & & 'Wrong working area.') DclSinQFT_B = r call sinqb(n, DclSinQFT_B, wk(idx)%array) end function end module !------------------------------------------------- ! FFT:偶数波数成分のみのcossine変換をおこなう.(6) !------------------------------------------------- module fftqcos use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitCosQFT(n,index) integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitCosQFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 3*n+15 allocate(wk(idx)%array(len)) call cosqi(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocCosQFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- function DclCosQFT_F(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclCosQFT_F if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosQFT_F', & & 'Wrong working area.') DclCosQFT_F = r call cosqf(n, DclCosQFT_F, wk(idx)%array) end function !--------------------------------------------------------- function DclCosQFT_B(r, index) real, intent(in), dimension(:) :: r integer, intent(in), optional :: index real, dimension(size(r)) :: DclCosQFT_B if(present(index)) then; idx = index else; idx = 1 ; end if n = size(r) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosQFT_B', & & 'Wrong working area.') DclCosQFT_B = r call cosqb(size(r), DclCosQFT_B, wk(idx)%array) end function end module !------------------------------------------------- ! 周期複素数データのフーリエ変換をおこなう (7) !------------------------------------------------- module fftcmplx use fft_work type(work), dimension(10), private :: wk contains !--------------------------------------------------------- subroutine DclInitComplexFFT(n,index) integer, intent(in) :: n !処理するデータの長さ integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitComplexFFT', & & 'Working area has been allocated already.') wk(idx)%n = n len = 4*n+15 allocate(wk(idx)%array(len)) call cffti(n, wk(idx)%array) end subroutine !--------------------------------------------------------- subroutine DclDeallocComplexFFT(index) integer, intent(in), optional :: index !作業領域番号 if(present(index)) then; idx = index else; idx = 1 ; end if deallocate(wk(idx)%array) end subroutine !--------------------------------------------------------- function DclComplexFFT_F(c,index) complex, intent(in), dimension(:) :: c integer, intent(in), optional :: index complex, dimension(size(c)) :: DclComplexFFT_F if(present(index)) then; idx = index else; idx = 1 ; end if n = size(c) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclComplexFFT_F', & & 'Wrong working area.') DclComplexFFT_F = c call cfftf(n, DclComplexFFT_F, wk(idx)%array) end function !--------------------------------------------------------- function DclComplexFFT_B(c,index) complex, intent(in), dimension(:) :: c integer, intent(in), optional :: index complex, dimension(size(c)) :: DclComplexFFT_B if(present(index)) then; idx = index else; idx = 1 ; end if n = size(c) if(n.ne.wk(idx)%n) call msgdmp('E', 'DclComplexFFT_B', & & 'Wrong working area.') DclComplexFFT_B = c call cfftb(n, DclComplexFFT_B, wk(idx)%array) end function end module !------------------------------------------------- ! fftlib Module !------------------------------------------------- module fftlib use fftreal use ffteasy use fftsin use fftcos use fftqsin use fftqcos use fftcmplx private :: work end module