Class | ffts |
In: |
ffts.f90
|
fft 関連のサブルーチン集
Subroutine : | |||
nx : | integer, intent(in)
| ||
a : | complex, intent(in), dimension(0:nx/2-1)
| ||
b : | real, intent(inout), dimension(0:nx-1)
| ||
prim : | character(1), optional, intent(in)
| ||
prim_fact : | integer, intent(in), optional, dimension(5)
| ||
omega_fix : | complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional
| ||
omegan_fix : | complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional
|
1 次元実 FFT 逆変換計算ルーチン
入力配列は複素数数で行い,実数配列を返す.
subroutine c2r_ffttp_1d( nx, a, b, prim, prim_fact, omega_fix, omegan_fix ) ! 1 次元実 FFT 逆変換計算ルーチン ! ! 入力配列は複素数数で行い,実数配列を返す. use Math_Const implicit none integer, intent(in) :: nx ! 1 次元データ要素数 complex, intent(in), dimension(0:nx/2-1) :: a ! 入力複素数データ配列 real, intent(inout), dimension(0:nx-1) :: b ! 出力実数データ配列 character(1), optional, intent(in) :: prim ! 素因数分解をするかどうか ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT. ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする. ! なお, 実 FFT 変換は現在, データ要素数が偶数の変換しか行わないように ! 実装されていることに注意 (奇数データで変換しようとするとエラーを返す). ! また, in 属性の引数は, データ数が半分になっていることに注意. ! 仕様として, 入力引数の n=0 虚部に nx/2 番目の実数データが入っているように ! データを渡す. ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する. integer, intent(in), optional, dimension(5) :: prim_fact ! prim = x のとき設定すると ! そのべきで FFT を行う. prim_fact=(/a,b,c,d/) : 2^a*3^b*5^c*7^d complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omega_fix ! 余素因数に対応する回転行列 complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omegan_fix ! 回転行列 complex, dimension(0:nx/2-1) :: c, d integer :: i integer, dimension(5) :: prim_num if(mod(nx,2)/=0)then write(*,*) "*** ERROR ***" write(*,*) "nx must be even number. Stop." stop end if !-- nx/2 個の複素数データ(独立なフーリエ係数)を同数の複素数データに置き換える. c(0)=a(0) do i=1,nx/2-1 c(i)=(a(i)+conjg(a(nx/2-i))) +(img*cos(2.0*pi*i/real(nx))-sin(2.0*pi*i/real(nx))) *(a(i)-conjg(a(nx/2-i))) end do !-- FFT 開始 if(present(prim))then if(present(prim_fact))then prim_num=prim_fact prim_num(1)=prim_num(1)-1 ! データの要素数が半分になっているため, ! 2 のべき数は 1 つ減っている. if(present(omegan_fix))then call ffttp_1d( nx/2, c, d, 'i', prim, prim_num, omega_fix(0:nx/2-1,0:nx/2-1), omegan_fix(0:nx/2-1,0:nx/2-1) ) else call ffttp_1d( nx/2, c, d, 'i', prim, prim_num ) end if else call ffttp_1d( nx/2, c, d, 'i', prim ) end if else call ffttp_1d( nx/2, c, d, 'i' ) end if !-- 変換後の配列を整理 do i=0,nx/2-1 b(2*i)=real(d(i)) b(2*i+1)=aimag(d(i)) end do end subroutine
Subroutine : | |||
nx : | integer, intent(in)
| ||
a : | complex, intent(in), dimension(0:nx-1)
| ||
b : | complex, intent(inout), dimension(0:nx-1)
| ||
csign : | character(1), intent(in)
| ||
prim : | character(1), optional, intent(in)
| ||
prim_fact : | integer, dimension(5), optional
| ||
omega_fix : | complex, dimension(0:nx-1,0:nx-1), intent(in), optional
| ||
omegan_fix : | complex, dimension(0:nx-1,0:nx-1), intent(in), optional
|
Temperton's FFT 1d の fft を csign をもとに, 正変換, 逆変換するルーチン
subroutine ffttp_1d( nx, a, b, csign, prim, prim_fact, omega_fix, omegan_fix ) ! Temperton's FFT ! ! 1d の fft を csign をもとに, 正変換, 逆変換するルーチン use Math_Const implicit none integer, intent(in) :: nx ! 入力配列の要素数 complex, intent(in), dimension(0:nx-1) :: a ! 入力配列 complex, intent(inout), dimension(0:nx-1) :: b ! 出力配列 character(1), intent(in) :: csign ! 正逆変換判定 [r=正変換, i=逆変換] character(1), optional, intent(in) :: prim ! 素因数分解をするかどうか ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT. ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする. ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する. integer, dimension(5), optional :: prim_fact ! prim = x のとき設定すると ! そのべきで FFT を行う. prim_fact=(/a,b,c,d,e/) : 2^a*3^b*5^c*7^d*e complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omega_fix ! 余素因数に対応する回転行列 complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omegan_fix ! 回転行列 integer, allocatable, dimension(:) :: l, m, n ! 要素等の作業用配列 complex :: ctmp integer :: stat, counter, base integer :: i, j, k, id, jd, kd ! do loop 用配列 integer, parameter, dimension(4) :: prim_dim=(/2, 3, 5, 7/) ! 素因数 integer, dimension(4) :: prim_num ! 各素因数のべき数 complex, allocatable, dimension(:,:) :: omega complex, dimension(0:nx-1,0:nx-1) :: omegan complex, dimension(0:nx-1) :: c ! tmp array complex, dimension(0:1,0:1) :: omega2 ! 波数 2 の回転行列 complex, dimension(0:2,0:2) :: omega3 ! 波数 3 の回転行列 complex, dimension(0:4,0:4) :: omega5 ! 波数 5 の回転行列 complex, dimension(0:6,0:6) :: omega7 ! 波数 7 の回転行列 base=nx prim_num=0 counter=0 do i=0,nx-1 b(i)=a(i) end do if(real(romega2(1,1))/=-1.0.or.aimag(romega2(1,1))/=1.0.or. real(iomega2(1,1))/=-1.0.or.aimag(iomega2(1,1))/=1.0)then ! math_const の rotate_array ルーチンが計算されていない場合, これを作動させる. call rotate_array() end if !-- 素因数分解する処理 if(present(prim))then if(prim=='o')then if(present(prim_fact))then do i=1,4 ! if(prim_fact(i)==0)then ! prim_num(i)=1 ! else prim_num(i)=prim_fact(i) ! end if counter=counter+prim_fact(i) end do base=prim_fact(5) else call prim_calc( nx, prim_num, base ) do i=1,4 counter=counter+prim_num(i) end do end if if(base==1)then counter=counter-1 end if if(counter/=0)then ! prim=='o' であっても, 素因数分解できなければ DFT に送る. allocate(l(counter+1)) allocate(m(counter+1)) allocate(n(counter+1)) stat=0 l=0 m=0 n=0 do i=1,4 if(prim_num(i)/=0)then select case(prim_dim(i)) case(2) n(stat+1:stat+prim_num(i))=2 case(3) n(stat+1:stat+prim_num(i))=3 case(5) n(stat+1:stat+prim_num(i))=5 case(7) n(stat+1:stat+prim_num(i))=7 end select stat=stat+prim_num(i) end if end do if(base/=1)then n(counter+1)=base end if end if end if end if !-- 回転行列を定義 if(counter/=0)then select case (csign) case ('r') do j=0,1 do i=0,1 omega2(i,j)=romega2(i,j) end do end do do j=0,2 do i=0,2 omega3(i,j)=romega3(i,j) end do end do do j=0,4 do i=0,4 omega5(i,j)=romega5(i,j) end do end do do j=0,6 do i=0,6 omega7(i,j)=romega7(i,j) end do end do case ('i') do j=0,1 do i=0,1 omega2(i,j)=iomega2(i,j) end do end do do j=0,2 do i=0,2 omega3(i,j)=iomega3(i,j) end do end do do j=0,4 do i=0,4 omega5(i,j)=iomega5(i,j) end do end do do j=0,6 do i=0,6 omega7(i,j)=iomega7(i,j) end do end do end select end if if(present(omegan_fix))then allocate(omega(0:base-1,0:base-1)) do j=0,base-1 do i=0,base-1 omega(i,j)=omega_fix(i,j) end do end do do j=0,nx-1 do i=0,nx-1 omegan(i,j)=omegan_fix(i,j) end do end do else allocate(omega(0:base-1,0:base-1)) call rotate_calc( nx, csign, (/prim_num(1), prim_num(2), prim_num(3), prim_num(4), base/), omega, omegan ) end if !-- FFT 開始 if(counter/=0)then !-- 係数行列定義 m(1)=1 l(1)=nx/(n(1)*m(1)) do i=2,counter+1 m(i)=m(i-1)*n(i-1) l(i)=nx/(n(i)*m(i)) end do !-- 変換行列 W の定義 do kd=1,counter+1 do jd=0,l(kd)-1 do id=0,n(kd)-1 do k=0,m(kd)-1 ctmp=b(jd*m(kd)+k) do j=1,n(kd)-1 select case(n(kd)) case(2) ctmp=ctmp+omega2(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k) case(3) ctmp=ctmp+omega3(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k) case(5) ctmp=ctmp+omega5(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k) case(7) ctmp=ctmp+omega7(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k) case default ctmp=ctmp+omega(id,j)*b(j*l(kd)*m(kd)+jd*m(kd)+k) end select end do c(jd*n(kd)*m(kd)+id*m(kd)+k)=ctmp*omegan(m(kd),(id*jd)) end do end do end do do id=0,nx-1 b(id)=c(id) end do end do else do j=0,nx-1 b(j)=a(0) do i=1,nx-1 b(j)=b(j)+a(i)*omegan(i,j) end do end do end if if(csign=='r')then do j=0,nx-1 b(j)=b(j)/real(nx) end do end if end subroutine
Subroutine : | |||
nx : | integer, intent(in)
| ||
ny : | integer, intent(in)
| ||
a : | complex, intent(in), dimension(0:nx-1,0:ny-1)
| ||
b : | complex, intent(inout), dimension(0:nx-1,0:ny-1)
| ||
csign : | character(1), intent(in)
| ||
prim : | character(1), optional, intent(in)
| ||
prim_fact : | integer, dimension(5), optional
| ||
omega_fix : | complex, dimension(0:nx-1,0:nx-1), intent(in), optional
| ||
omegan_fix : | complex, dimension(0:nx-1,0:nx-1), intent(in), optional
|
Temperton's FFT (2d ver) 2d の fft を csign をもとに, 正変換, 逆変換するルーチン
subroutine ffttp_2d( nx, ny, a, b, csign, prim, prim_fact, omega_fix, omegan_fix ) ! Temperton's FFT (2d ver) ! ! 2d の fft を csign をもとに, 正変換, 逆変換するルーチン use Math_Const implicit none integer, intent(in) :: nx ! 入力配列の要素数 1 integer, intent(in) :: ny ! 入力配列の要素数 2 complex, intent(in), dimension(0:nx-1,0:ny-1) :: a ! 入力配列 complex, intent(inout), dimension(0:nx-1,0:ny-1) :: b ! 出力配列 character(1), intent(in) :: csign ! 正逆変換判定 [r=正変換, i=逆変換] character(1), optional, intent(in) :: prim ! 素因数分解をするかどうか ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT. ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする. ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する. integer, dimension(5), optional :: prim_fact ! prim = x のとき設定すると ! そのべきで FFT を行う. prim_fact=(/a,b,c,d/) : 2^a*3^b*5^c*7^d ! ここで, prim_fact で設定するべき数は nx 方向のべき数であることに注意. complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omega_fix ! 余素因数に対応する回転行列 complex, dimension(0:nx-1,0:nx-1), intent(in), optional :: omegan_fix ! 回転行列 integer :: i if(present(prim))then if(present(prim_fact))then if(present(omegan_fix))then do i=0,ny-1 call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign, prim, prim_fact, omega_fix(0:nx-1,0:nx-1), omegan_fix(0:nx-1,0:nx-1) ) end do else do i=0,ny-1 call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign, prim, prim_fact ) end do end if else do i=0,ny-1 call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign, prim ) end do end if else do i=0,ny-1 call ffttp_1d( nx, a(0:nx-1,i), b(0:nx-1,i), csign ) end do end if end subroutine
Subroutine : | |||
n : | integer, intent(in)
| ||
factor : | integer, intent(inout), dimension(4)
| ||
resid : | integer, intent(inout)
|
2,3,5,7 についての素因数分解を行う.
subroutine prim_calc( n, factor, resid ) ! 2,3,5,7 についての素因数分解を行う. implicit none integer, intent(in) :: n ! 分解する数値 integer, intent(inout), dimension(4) :: factor ! 分解したときの各素因数 ! このとき, factor(1:4)=(/a,b,c,d/) であり, factor(5)=e とすると, ! n=2^a*3^b*5^c*7^d*e という式になっている. integer, intent(inout) :: resid ! residual prim factor integer :: i, base integer, dimension(4) :: prim_dim prim_dim=(/2,3,5,7/) base=n factor=(/0,0,0,0/) do i=1,4 do while(mod(base,prim_dim(i))==0) base=base/prim_dim(i) factor(i)=factor(i)+1 end do end do resid=base end subroutine
Subroutine : | |||
nx : | integer, intent(in)
| ||
a : | real, intent(in), dimension(0:nx-1)
| ||
b : | complex, intent(inout), dimension(0:nx/2-1)
| ||
prim : | character(1), optional, intent(in)
| ||
prim_fact : | integer, dimension(5), intent(in), optional
| ||
omega_fix : | complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional
| ||
omegan_fix : | complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional
|
1 次元実 FFT 計算ルーチン
入力配列は実数で行い, 複素数配列を返す.
subroutine r2c_ffttp_1d( nx, a, b, prim, prim_fact, omega_fix, omegan_fix ) ! 1 次元実 FFT 計算ルーチン ! ! 入力配列は実数で行い, 複素数配列を返す. use Math_Const implicit none integer, intent(in) :: nx ! 1 次元データ要素数 real, intent(in), dimension(0:nx-1) :: a ! 入力実数データ配列 complex, intent(inout), dimension(0:nx/2-1) :: b ! 出力複素数データ配列 character(1), optional, intent(in) :: prim ! 素因数分解をするかどうか ! [o=分解する, x=分解しない] default=分解しない. その場合は, 通常の DFT. ! 素因数分解する場合, nx=2^a*3^b*5^c*7^d までしか分解しないようにする. ! なお, 実 FFT 変換は現在, データ要素数が偶数の変換しか行わないように ! 実装されていることに注意 (奇数データで変換しようとするとエラーを返す). ! また, inout 属性の引数は, データ数が半分になっていることに注意. ! 仕様として, nx/2 番目の実数データは, n=0 番目の虚部に格納されて ! 出力配列に渡されることに注意. ! ただし, 'x' のときでも, prim_fact が設定されていれば, そのべきで FFT する. integer, dimension(5), intent(in), optional :: prim_fact ! prim = x のとき設定すると ! そのべきで FFT を行う. prim_fact=(/a,b,c,d/) : 2^a*3^b*5^c*7^d ! ここで入力するべき数は nx のべき数である. complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omega_fix ! 余素因数に対応する回転行列 complex, dimension(0:nx/2-1,0:nx/2-1), intent(in), optional :: omegan_fix ! 回転行列 complex, dimension(0:nx/2-1) :: c, d integer :: i integer, dimension(5) :: prim_num if(mod(nx,2)/=0)then write(*,*) "*** ERROR ***" write(*,*) "nx must be even number. Stop." stop end if !-- nx 個の実数データを nx/2 個の複素数データに置き換える. do i=0,nx/2-1 c(i)=a(2*i)+img*a(2*i+1) end do !-- FFT 開始 if(present(prim))then if(present(prim_fact))then prim_num=prim_fact prim_num(1)=prim_num(1)-1 ! データの要素数が半分になっているため, ! 2 のべき数は 1 つ減っている. if(present(omegan_fix))then call ffttp_1d( nx/2, c, d, 'r', prim, prim_num, omega_fix(0:nx/2-1,0:nx/2-1), omegan_fix(0:nx/2-1,0:nx/2-1) ) else call ffttp_1d( nx/2, c, d, 'r', prim, prim_num ) end if else call ffttp_1d( nx/2, c, d, 'r', prim ) end if else call ffttp_1d( nx/2, c, d, 'r' ) end if !-- 変換後の配列を整理 !-- b(N) の実部は b(0) の虚部に組み込むことにする. !-- b(k) の計算で 0.25 をかけるのは, 上の fft で 2/N で規格化しており, !-- もとの計算では, 1/2N で規格化しなければならないので, !-- 1/4 をかけることで, 2/N -> 1/2N で規格化したことになる. !-- b(0) にかかっている係数もその類. b(0)=0.5*((real(d(0))+aimag(d(0)))+img*(real(d(0))-aimag(d(0)))) do i=1,nx/2-1 b(i)=0.25*((conjg(d(nx/2-i))+d(i)) -(sin(2.0*pi*real(i)/real(nx))+img*cos(2.0*pi*real(i)/real(nx))) *(d(i)-conjg(d(nx/2-i)))) end do end subroutine
Subroutine : | |||
nx : | integer, intent(in)
| ||
csign : | character(1), intent(in)
| ||
prim_fact : | integer, dimension(5)
| ||
omega : | complex, dimension(0:prim_fact(5)-1,0:prim_fact(5)-1),
intent(inout)
| ||
omegan : | complex, dimension(0:nx-1,0:nx-1), intent(inout)
|
FFT に使用する回転行列を計算する.
subroutine rotate_calc( nx, csign, prim_fact, omega, omegan ) ! FFT に使用する回転行列を計算する. use Math_Const implicit none integer, intent(in) :: nx ! 入力配列の要素数 character(1), intent(in) :: csign ! 正逆変換判定 [r=正変換, i=逆変換] integer, dimension(5) :: prim_fact ! prim = x のとき設定すると ! そのべきで FFT を行う. prim_fact=(/a,b,c,d,e/) : 2^a*3^b*5^c*7^d*e complex, dimension(0:prim_fact(5)-1,0:prim_fact(5)-1), intent(inout) :: omega ! 余因数に対応する回転行列 complex, dimension(0:nx-1,0:nx-1), intent(inout) :: omegan ! 回転行列 integer :: counter, base integer :: i, j ! do loop 用配列 base=prim_fact(5) counter=0 do i=1,4 counter=counter+prim_fact(i) end do !-- 回転行列を定義 select case(csign) case('r') if(counter/=0)then if(base/=1)then do j=0,base-1 do i=0,base-1 omega(i,j)=cos(2.0*pi*i*j/real(base))-img*sin(2.0*pi*i*j/real(base)) end do end do end if do j=0,nx-1 do i=0,nx-1 omegan(i,j)=cos(2.0*pi*i*j/real(nx))-img*sin(2.0*pi*i*j/real(nx)) end do end do else do j=0,nx-1 do i=0,nx-1 omega(i,j)=cos(2.0*pi*i*j/real(nx))-img*sin(2.0*pi*i*j/real(nx)) omegan(i,j)=omega(i,j) end do end do end if case('i') if(counter/=0)then if(base/=1)then do j=0,base-1 do i=0,base-1 omega(i,j)=cos(2.0*pi*i*j/real(base))+img*sin(2.0*pi*i*j/real(base)) end do end do end if do j=0,nx-1 do i=0,nx-1 omegan(i,j)=cos(2.0*pi*i*j/real(nx))+img*sin(2.0*pi*i*j/real(nx)) end do end do else do j=0,nx-1 do i=0,nx-1 omega(i,j)=cos(2.0*pi*i*j/real(nx))+img*sin(2.0*pi*i*j/real(nx)) omegan(i,j)=omega(i,j) end do end do end if case default write(*,*) "******** ERROR : csign is bad. **********" write(*,*) "Stop!" stop end select end subroutine