fftsg3d.f

2维fft程序

! Fast Fourier/Cosine/Sine Transform
!     dimension   :three
!     data length :power of 2
!     decimation  :frequency
!     radix       :split-radix, row-column
!     data        :inplace
!     table       :use
! subroutines
!     cdft3d: Complex Discrete Fourier Transform
!     rdft3d: Real Discrete Fourier Transform
!     ddct3d: Discrete Cosine Transform
!     ddst3d: Discrete Sine Transform
! necessary package
!     fftsg.f  : 1D-FFT package
!
!
! -------- Complex DFT (Discrete Fourier Transform) --------
!     [definition]
!         <case1>
!             X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               x(j1,j2,j3) * 
!                               exp(2*pi*i*j1*k1/n1) * 
!                               exp(2*pi*i*j2*k2/n2) * 
!                               exp(2*pi*i*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2>
!             X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               x(j1,j2,j3) * 
!                               exp(-2*pi*i*j1*k1/n1) * 
!                               exp(-2*pi*i*j2*k2/n2) * 
!                               exp(-2*pi*i*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         2*n1   :data length (integer)
!                 n1 >= 1, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 1, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 1, n3 = power of 2
!         a(0:2*n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 input data
!                     a(2*j1,j2,j3) = Re(x(j1,j2,j3)), 
!                     a(2*j1+1,j2,j3) = Im(x(j1,j2,j3)), 
!                     0<=j1<n1, 0<=j2<n2, 0<=j3<n3
!                 output data
!                     a(2*k1,k2,k3) = Re(X(k1,k2,k3)), 
!                     a(2*k1+1,k2,k3) = Im(X(k1,k2,k3)), 
!                     0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         t(0:*) :work area (real*8)
!                 length of t >= max(8*n2, 8*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1, n2, n3))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1/2, n2/2, n3/2)
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w)
!         is 
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, 2 * n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (1.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
! -------- Real DFT / Inverse of Real DFT --------
!     [definition]
!         <case1> RDFT
!             R(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!             I(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2> IRDFT (excluding scale)
!             a(k1,k2,k3) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               (R(j1,j2,j3) * 
!                               cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3) + 
!                               I(j1,j2,j3) * 
!                               sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3)), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         (notes: R(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = R(k1,k2,k3), 
!                 I(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = -I(k1,k2,k3), 
!                 0<=k1<n1, 0<=k2<n2, 0<=k3<n3)
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         n1     :data length (integer)
!                 n1 >= 2, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 2, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 2, n3 = power of 2
!         a(0:n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 <case1>
!                     output data
!                         a(2*k1,k2,k3) = R(k1,k2,k3)
!                                 = R(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), 
!                         a(2*k1+1,k2,k3) = I(k1,k2,k3)
!                                 = -I(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), 
!                             0<k1<n1/2, 0<=k2<n2, 0<=k3<n3, 
!                         a(0,k2,k3) = R(0,k2,k3)
!                                    = R(0,n2-k2,mod(n3-k3,n3)), 
!                         a(1,k2,k3) = I(0,k2,k3)
!                                    = -I(0,n2-k2,mod(n3-k3,n3)), 
!                         a(1,n2-k2,k3) = R(n1/2,k2,mod(n3-k3,n3))
!                                       = R(n1/2,n2-k2,k3), 
!                         a(0,n2-k2,k3) = -I(n1/2,k2,mod(n3-k3,n3))
!                                       = I(n1/2,n2-k2,k3), 
!                             0<k2<n2/2, 0<=k3<n3, 
!                         a(0,0,k3) = R(0,0,k3)
!                                   = R(0,0,n3-k3), 
!                         a(1,0,k3) = I(0,0,k3)
!                                   = -I(0,0,n3-k3), 
!                         a(0,n2/2,k3) = R(0,n2/2,k3)
!                                      = R(0,n2/2,n3-k3), 
!                         a(1,n2/2,k3) = I(0,n2/2,k3)
!                                      = -I(0,n2/2,n3-k3), 
!                         a(1,0,n3-k3) = R(n1/2,0,k3)
!                                      = R(n1/2,0,n3-k3), 
!                         a(0,0,n3-k3) = -I(n1/2,0,k3)
!                                      = I(n1/2,0,n3-k3), 
!                         a(1,n2/2,n3-k3) = R(n1/2,n2/2,k3)
!                                         = R(n1/2,n2/2,n3-k3), 
!                         a(0,n2/2,n3-k3) = -I(n1/2,n2/2,k3)
!                                         = I(n1/2,n2/2,n3-k3), 
!                             0<k3<n3/2, 
!                         a(0,0,0) = R(0,0,0), 
!                         a(1,0,0) = R(n1/2,0,0), 
!                         a(0,0,n3/2) = R(0,0,n3/2), 
!                         a(1,0,n3/2) = R(n1/2,0,n3/2), 
!                         a(0,n2/2,0) = R(0,n2/2,0), 
!                         a(1,n2/2,0) = R(n1/2,n2/2,0), 
!                         a(0,n2/2,n3/2) = R(0,n2/2,n3/2), 
!                         a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2)
!                 <case2>
!                     input data
!                         a(2*j1,j2,j3) = R(j1,j2,j3)
!                                 = R(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), 
!                         a(2*j1+1,j2,j3) = I(j1,j2,j3)
!                                 = -I(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), 
!                             0<j1<n1/2, 0<=j2<n2, 0<=j3<n3, 
!                         a(0,j2,j3) = R(0,j2,j3)
!                                    = R(0,n2-j2,mod(n3-j3,n3)), 
!                         a(1,j2,j3) = I(0,j2,j3)
!                                    = -I(0,n2-j2,mod(n3-j3,n3)), 
!                         a(1,n2-j2,j3) = R(n1/2,j2,mod(n3-j3,n3))
!                                       = R(n1/2,n2-j2,j3), 
!                         a(0,n2-j2,j3) = -I(n1/2,j2,mod(n3-j3,n3))
!                                       = I(n1/2,n2-j2,j3), 
!                             0<j2<n2/2, 0<=j3<n3, 
!                         a(0,0,j3) = R(0,0,j3)
!                                   = R(0,0,n3-j3), 
!                         a(1,0,j3) = I(0,0,j3)
!                                   = -I(0,0,n3-j3), 
!                         a(0,n2/2,j3) = R(0,n2/2,j3)
!                                      = R(0,n2/2,n3-j3), 
!                         a(1,n2/2,j3) = I(0,n2/2,j3)
!                                      = -I(0,n2/2,n3-j3), 
!                         a(1,0,n3-j3) = R(n1/2,0,j3)
!                                      = R(n1/2,0,n3-j3), 
!                         a(0,0,n3-j3) = -I(n1/2,0,j3)
!                                      = I(n1/2,0,n3-j3), 
!                         a(1,n2/2,n3-j3) = R(n1/2,n2/2,j3)
!                                         = R(n1/2,n2/2,n3-j3), 
!                         a(0,n2/2,n3-j3) = -I(n1/2,n2/2,j3)
!                                         = I(n1/2,n2/2,n3-j3), 
!                             0<j3<n3/2, 
!                         a(0,0,0) = R(0,0,0), 
!                         a(1,0,0) = R(n1/2,0,0), 
!                         a(0,0,n3/2) = R(0,0,n3/2), 
!                         a(1,0,n3/2) = R(n1/2,0,n3/2), 
!                         a(0,n2/2,0) = R(0,n2/2,0), 
!                         a(1,n2/2,0) = R(n1/2,n2/2,0), 
!                         a(0,n2/2,n3/2) = R(0,n2/2,n3/2), 
!                         a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2)
!                 ---- output ordering ----
!                     call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!                     call rdft3dsort(n1max, n2max, n1, n2, n3, 1, a)
!                     ! stored data is a(0:n1-1,0:n2-1,0:n3+1):
!                     ! a(2*k1,k2,k3) = R(k1,k2,k3), 
!                     ! a(2*k1+1,k2,k3) = I(k1,k2,k3), 
!                     ! 0<=k1<=n1/2, 0<=k2<n2, 0<=k3<n3.
!                     ! the stored data is larger than the input data!
!                 ---- input ordering ----
!                     call rdft3dsort(n1max, n2max, n1, n2, n3, -1, a)
!                     call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!         t(0:*) :work area (real*8)
!                 length of t >= max(8*n2, 8*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1/2, n2, n3))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1/4, n2/2, n3/2) + n1/4
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         is 
!             call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (2.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
!     [definition]
!         <case1> IDCT (excluding scale)
!             C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               cos(pi*j1*(k1+1/2)/n1) * 
!                               cos(pi*j2*(k2+1/2)/n2) * 
!                               cos(pi*j3*(k3+1/2)/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2> DCT
!             C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               cos(pi*(j1+1/2)*k1/n1) * 
!                               cos(pi*(j2+1/2)*k2/n2) * 
!                               cos(pi*(j3+1/2)*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         n1     :data length (integer)
!                 n1 >= 2, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 2, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 2, n3 = power of 2
!         a(0:n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 output data
!                     a(k1,k2,k3) = C(k1,k2,k3), 
!                         0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         t(0:*) :work area (real*8)
!                 length of t >= max(4*n2, 4*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1/2, n2/2, n3/2))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1*3/2, n2*3/2, n3*3/2)
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!         is 
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     a(0, j2, j3) = a(0, j2, j3) * 0.5d0
!                 end do
!                 do j1 = 0, n1 - 1
!                     a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0
!                 end do
!             end do
!             do j2 = 0, n2 - 1
!                 do j1 = 0, n1 - 1
!                     a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0
!                 end do
!             end do
!             call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1