| ! Fast Fourier/Cosine/Sine Transform |
| ! dimension :two |
| ! data length :power of 2 |
| ! decimation :frequency |
| ! radix :4, 2, row-column |
| ! data :inplace |
| ! table :use |
| ! subroutines |
| ! cdft2d: Complex Discrete Fourier Transform |
| ! rdft2d: Real Discrete Fourier Transform |
| ! ddct2d: Discrete Cosine Transform |
| ! ddst2d: Discrete Sine Transform |
| ! |
| ! |
| ! -------- Complex DFT (Discrete Fourier Transform) -------- |
| ! [definition] |
| ! <case1> |
| ! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) * |
| ! exp(2*pi*i*j1*k1/n1) * |
| ! exp(2*pi*i*j2*k2/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! <case2> |
| ! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) * |
| ! exp(-2*pi*i*j1*k1/n1) * |
| ! exp(-2*pi*i*j2*k2/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! (notes: sum_j=0^n-1 is a summation from j=0 to n-1) |
| ! [usage] |
| ! <case1> |
| ! ip(0) = 0 ! first time only |
| ! call cdft2d(n1max, 2*n1, n2, 1, a, ip, w) |
| ! <case2> |
| ! ip(0) = 0 ! first time only |
| ! call cdft2d(n1max, 2*n1, n2, -1, a, ip, w) |
| ! [parameters] |
| ! n1max :row size of the 2D array (integer) |
| ! 2*n1 :data length (integer) |
| ! n1 >= 1, n1 = power of 2 |
| ! n2 :data length (integer) |
| ! n2 >= 1, n2 = power of 2 |
| ! a(0:2*n1-1,0:n2-1) |
| ! :input/output data (real*8) |
| ! input data |
| ! a(2*j1,j2) = Re(x(j1,j2)), |
| ! a(2*j1+1,j2) = Im(x(j1,j2)), |
| ! 0<=j1<n1, 0<=j2<n2 |
| ! output data |
| ! a(2*k1,k2) = Re(X(k1,k2)), |
| ! a(2*k1+1,k2) = Im(X(k1,k2)), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! ip(0:*):work area for bit reversal (integer) |
| ! length of ip >= 2+sqrt(n) |
| ! (n = max(n1, n2)) |
| ! 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) |
| ! w(),ip() are initialized if ip(0) = 0. |
| ! [remark] |
| ! Inverse of |
| ! call cdft2d(n1max, 2*n1, n2, -1, a, ip, w) |
| ! is |
| ! call cdft2d(n1max, 2*n1, n2, 1, a, ip, w) |
| ! do j2 = 0, n2 - 1 |
| ! do j1 = 0, 2 * n1 - 1 |
| ! a(j1, j2) = a(j1, j2) * (1.0d0 / (n1 * n2)) |
| ! end do |
| ! end do |
| ! . |
| ! |
| ! |
| ! -------- Real DFT / Inverse of Real DFT -------- |
| ! [definition] |
| ! <case1> RDFT |
| ! R(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * |
| ! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! I(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * |
| ! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! <case2> IRDFT (excluding scale) |
| ! a(k1,k2) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 |
| ! (R(j1,j2) * |
| ! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2) + |
| ! I(j1,j2) * |
| ! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2)), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! (notes: R(n1-k1,n2-k2) = R(k1,k2), |
| ! I(n1-k1,n2-k2) = -I(k1,k2), |
| ! R(n1-k1,0) = R(k1,0), |
| ! I(n1-k1,0) = -I(k1,0), |
| ! R(0,n2-k2) = R(0,k2), |
| ! I(0,n2-k2) = -I(0,k2), |
| ! 0<k1<n1, 0<k2<n2) |
| ! [usage] |
| ! <case1> |
| ! ip(0) = 0 ! first time only |
| ! call rdft2d(n1max, n1, n2, 1, a, ip, w) |
| ! <case2> |
| ! ip(0) = 0 ! first time only |
| ! call rdft2d(n1max, n1, n2, -1, a, ip, w) |
| ! [parameters] |
| ! n1max :row size of the 2D array (integer) |
| ! n1 :data length (integer) |
| ! n1 >= 2, n1 = power of 2 |
| ! n2 :data length (integer) |
| ! n2 >= 2, n2 = power of 2 |
| ! a(0:n1-1,0:n2-1) |
| ! :input/output data (real*8) |
| ! <case1> |
| ! output data |
| ! a(2*k1,k2) = R(k1,k2) = R(n1-k1,n2-k2), |
| ! a(2*k1+1,k2) = I(k1,k2) = -I(n1-k1,n2-k2), |
| ! 0<k1<n1/2, 0<k2<n2, |
| ! a(2*k1,0) = R(k1,0) = R(n1-k1,0), |
| ! a(2*k1+1,0) = I(k1,0) = -I(n1-k1,0), |
| ! 0<k1<n1/2, |
| ! a(0,k2) = R(0,k2) = R(0,n2-k2), |
| ! a(1,k2) = I(0,k2) = -I(0,n2-k2), |
| ! a(1,n2-k2) = R(n1/2,k2) = R(n1/2,n2-k2), |
| ! a(0,n2-k2) = -I(n1/2,k2) = I(n1/2,n2-k2), |
| ! 0<k2<n2/2, |
| ! a(0,0) = R(0,0), |
| ! a(1,0) = R(n1/2,0), |
| ! a(0,n2/2) = R(0,n2/2), |
| ! a(1,n2/2) = R(n1/2,n2/2) |
| ! <case2> |
| ! input data |
| ! a(2*j1,j2) = R(j1,j2) = R(n1-j1,n2-j2), |
| ! a(2*j1+1,j2) = I(j1,j2) = -I(n1-j1,n2-j2), |
| ! 0<j1<n1/2, 0<j2<n2, |
| ! a(2*j1,0) = R(j1,0) = R(n1-j1,0), |
| ! a(2*j1+1,0) = I(j1,0) = -I(n1-j1,0), |
| ! 0<j1<n1/2, |
| ! a(0,j2) = R(0,j2) = R(0,n2-j2), |
| ! a(1,j2) = I(0,j2) = -I(0,n2-j2), |
| ! a(1,n2-j2) = R(n1/2,j2) = R(n1/2,n2-j2), |
| ! a(0,n2-j2) = -I(n1/2,j2) = I(n1/2,n2-j2), |
| ! 0<j2<n2/2, |
| ! a(0,0) = R(0,0), |
| ! a(1,0) = R(n1/2,0), |
| ! a(0,n2/2) = R(0,n2/2), |
| ! a(1,n2/2) = R(n1/2,n2/2) |
| ! ip(0:*):work area for bit reversal (integer) |
| ! length of ip >= 2+sqrt(n) |
| ! (n = max(n1/2, n2)) |
| ! 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) + n1/4 |
| ! w(),ip() are initialized if ip(0) = 0. |
| ! [remark] |
| ! Inverse of |
| ! call rdft2d(n1max, n1, n2, 1, a, ip, w) |
| ! is |
| ! call rdft2d(n1max, n1, n2, -1, a, ip, w) |
| ! do j2 = 0, n2 - 1 |
| ! do j1 = 0, n1 - 1 |
| ! a(j1, j2) = a(j1, j2) * (2.0d0 / (n1 * n2)) |
| ! end do |
| ! end do |
| ! . |
| ! |
| ! |
| ! -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- |
| ! [definition] |
| ! <case1> IDCT (excluding scale) |
| ! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * |
| ! cos(pi*j1*(k1+1/2)/n1) * |
| ! cos(pi*j2*(k2+1/2)/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! <case2> DCT |
| ! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * |
| ! cos(pi*(j1+1/2)*k1/n1) * |
| ! cos(pi*(j2+1/2)*k2/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! [usage] |
| ! <case1> |
| ! ip(0) = 0 ! first time only |
| ! call ddct2d(n1max, n1, n2, 1, a, t, ip, w) |
| ! <case2> |
| ! ip(0) = 0 ! first time only |
| ! call ddct2d(n1max, n1, n2, -1, a, t, ip, w) |
| ! [parameters] |
| ! n1max :row size of the 2D array (integer) |
| ! n1 :data length (integer) |
| ! n1 >= 2, n1 = power of 2 |
| ! n2 :data length (integer) |
| ! n2 >= 2, n2 = power of 2 |
| ! a(0:n1-1,0:n2-1) |
| ! :input/output data (real*8) |
| ! output data |
| ! a(k1,k2) = C(k1,k2), 0<=k1<n1, 0<=k2<n2 |
| ! t(0:n1-1,0:n2-1) |
| ! :work area (real*8) |
| ! ip(0:*):work area for bit reversal (integer) |
| ! length of ip >= 2+sqrt(n) |
| ! (n = max(n1/2, n2)) |
| ! 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) + max(n1, n2) |
| ! w(),ip() are initialized if ip(0) = 0. |
| ! [remark] |
| ! Inverse of |
| ! call ddct2d(n1max, n1, n2, -1, a, t, ip, w) |
| ! is |
| ! do j1 = 0, n1 - 1 |
| ! a(j1, 0) = a(j1, 0) * 0.5d0 |
| ! end do |
| ! do j2 = 0, n2 - 1 |
| ! a(0, j2) = a(0, j2) * 0.5d0 |
| ! end do |
| ! call ddct2d(n1max, n1, n2, 1, a, t, ip, w) |
| ! do j2 = 0, n2 - 1 |
| ! do j1 = 0, n1 - 1 |
| ! a(j1, j2) = a(j1, j2) * (4.0d0 / (n1 * n2)) |
| ! end do |
| ! end do |
| ! . |
| ! |
| ! |
| ! -------- DST (Discrete Sine Transform) / Inverse of DST -------- |
| ! [definition] |
| ! <case1> IDST (excluding scale) |
| ! S(k1,k2) = sum_j1=1^n1 sum_j2=1^n2 A(j1,j2) * |
| ! sin(pi*j1*(k1+1/2)/n1) * |
| ! sin(pi*j2*(k2+1/2)/n2), |
| ! 0<=k1<n1, 0<=k2<n2 |
| ! <case2> DST |
| ! S(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * |
| ! sin(pi*(j1+1/2)*k1/n1) * |
| ! sin(pi*(j2+1/2)*k2/n2), |
| ! 0<k1<=n1, 0<k2<=n2 |
| ! [usage] |
| ! <case1> |
| ! ip(0) = 0 ! first time only |
| ! call ddst2d(n1max, n1, n2, 1, a, t, ip, w) |
| ! <case2> |
| ! ip(0) = 0 ! first time only |
| ! call ddst2d(n1max, n1, n2, -1, a, t, ip, w) |
| ! [parameters] |
| ! n1max :row size of the 2D array (integer) |
| ! n1 :data length (integer) |
| ! n1 >= 2, n1 = power of 2 |
| ! n2 :data length (integer) |
| ! n2 >= 2, n2 = power of 2 |
| ! a(0:n1-1,0:n2-1) |
| ! :input/output data (real*8) |
| ! <case1> |
| ! input data |
| ! a(j1,j2) = A(j1,j2), 0<j1<n1, 0<j2<n2, |
| ! a(j1,0) = A(j1,n2), 0<j1<n1, |
| ! a(0,j2) = A(n1,j2), 0<j2<n2, |
| ! a(0,0) = A(n1,n2) |
| ! (i.e. A(j1,j2) = a(mod(j1,n1),mod(j2,n2))) |
| ! output data |
| ! a(k1,k2) = S(k1,k2), 0<=k1<n1, 0<=k2<n2 |
| ! <case2> |
| ! output data |
| ! a(k1,k2) = S(k1,k2), 0<k1<n1, 0<k2<n2, |
| ! a(k1,0) = S(k1,n2), 0<k1<n1, |
| ! a(0,k2) = S(n1,k2), 0<k2<n2, |
| ! a(0,0) = S(n1,n2) |
| ! (i.e. S(k1,k2) = a(mod(k1,n1),mod(k2,n2))) |
| ! t(0:n1-1,0:n2-1) |
| ! :work area (real*8) |
| ! ip(0:*):work area for bit reversal (integer) |
| ! length of ip >= 2+sqrt(n) |
| ! (n = max(n1/2, n2)) |
| ! 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) + max(n1, n2) |
| ! w(),ip() are initialized if ip(0) = 0. |
| ! [remark] |
| ! Inverse of |
| ! call ddst2d(n1max, n1, n2, -1, a, t, ip, w) |
| ! is |
| ! do j1 = 0, n1 - 1 |
| ! a(j1, 0) = a(j1, 0) * 0.5d0 |
| ! end do |
| ! do j2 = 0, n2 - 1 |
| ! a(0, j2) = a(0, j2) * 0.5d0 |
| ! end do |
| ! call ddst2d(n1max, n1, n2, 1, a, t, ip, w) |
| ! do j2 = 0, n2 - 1 |
| ! do j1 = 0, n1 - 1 |
| ! a(j1, j2) = a(j1, j2) * (4.0d0 / (n1 * n2)) |
| ! end do |
| ! end do |
| ! . |
| ! |
| ! |
| subroutine cdft2d(n1max, n1, n2, isgn, a, ip, w) |
| integer n1max, n1, n2, isgn, ip(0 : *), n |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *) |
| n = max(n1, 2 * n2) |
| if (n .gt. 4 * ip(0)) then |
| call makewt(n / 4, ip, w) |
| end if |
| if (n1 .gt. 4) then |
| call bitrv2row(n1max, n1, n2, ip(2), a) |
| end if |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), a) |
| end if |
| if (isgn .lt. 0) then |
| call cftfrow(n1max, n1, n2, a, w) |
| call cftfcol(n1max, n1, n2, a, w) |
| else |
| call cftbrow(n1max, n1, n2, a, w) |
| call cftbcol(n1max, n1, n2, a, w) |
| end if |
| end |
| ! |
| subroutine rdft2d(n1max, n1, n2, isgn, a, ip, w) |
| integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, n2h, i, j |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *), xi |
| n = max(n1, 2 * n2) |
| nw = ip(0) |
| if (n .gt. 4 * nw) then |
| nw = n / 4 |
| call makewt(nw, ip, w) |
| end if |
| nc = ip(1) |
| if (n1 .gt. 4 * nc) then |
| nc = n1 / 4 |
| call makect(nc, ip, w(nw)) |
| end if |
| n2h = n2 / 2 |
| if (isgn .lt. 0) then |
| do i = 1, n2h - 1 |
| j = n2 - i |
| xi = a(0, i) - a(0, j) |
| a(0, i) = a(0, i) + a(0, j) |
| a(0, j) = xi |
| xi = a(1, j) - a(1, i) |
| a(1, i) = a(1, i) + a(1, j) |
| a(1, j) = xi |
| end do |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), a) |
| end if |
| call cftfcol(n1max, n1, n2, a, w) |
| do i = 0, n2 - 1 |
| a(1, i) = 0.5d0 * (a(0, i) - a(1, i)) |
| a(0, i) = a(0, i) - a(1, i) |
| end do |
| if (n1 .gt. 4) then |
| call rftfrow(n1max, n1, n2, a, nc, w(nw)) |
| call bitrv2row(n1max, n1, n2, ip(2), a) |
| end if |
| call cftfrow(n1max, n1, n2, a, w) |
| else |
| if (n1 .gt. 4) then |
| call bitrv2row(n1max, n1, n2, ip(2), a) |
| end if |
| call cftbrow(n1max, n1, n2, a, w) |
| if (n1 .gt. 4) then |
| call rftbrow(n1max, n1, n2, a, nc, w(nw)) |
| end if |
| do i = 0, n2 - 1 |
| xi = a(0, i) - a(1, i) |
| a(0, i) = a(0, i) + a(1, i) |
| a(1, i) = xi |
| end do |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), a) |
| end if |
| call cftbcol(n1max, n1, n2, a, w) |
| do i = 1, n2h - 1 |
| j = n2 - i |
| a(0, j) = 0.5d0 * (a(0, i) - a(0, j)) |
| a(0, i) = a(0, i) - a(0, j) |
| a(1, j) = 0.5d0 * (a(1, i) + a(1, j)) |
| a(1, i) = a(1, i) - a(1, j) |
| end do |
| end if |
| end |
| ! |
| subroutine ddct2d(n1max, n1, n2, isgn, a, t, ip, w) |
| integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, n1h, n2h, |
| & i, ix, ic, j, jx, jc |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), |
| & t(0 : n1max - 1, 0 : n2 - 1), w(0 : *), xi |
| n = max(n1, 2 * n2) |
| nw = ip(0) |
| if (n .gt. 4 * nw) then |
| nw = n / 4 |
| call makewt(nw, ip, w) |
| end if |
| nc = ip(1) |
| if (n1 .gt. nc .or. n2 .gt. nc) then |
| nc = max(n1, n2) |
| call makect(nc, ip, w(nw)) |
| end if |
| n1h = n1 / 2 |
| n2h = n2 / 2 |
| if (isgn .ge. 0) then |
| do i = 0, n2 - 1 |
| do j = 1, n1h - 1 |
| jx = 2 * j |
| t(jx, i) = a(j, i) |
| t(jx + 1, i) = a(n1 - j, i) |
| end do |
| end do |
| t(0, 0) = a(0, 0) |
| t(1, 0) = a(n1h, 0) |
| t(0, n2h) = a(0, n2h) |
| t(1, n2h) = a(n1h, n2h) |
| do i = 1, n2h - 1 |
| ic = n2 - i |
| t(0, i) = a(0, i) |
| t(1, ic) = a(n1h, i) |
| t(1, i) = a(0, ic) |
| t(0, ic) = a(n1h, ic) |
| end do |
| call dctfsub(n1max, n1, n2, t, nc, w(nw)) |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), t) |
| end if |
| call cftfcol(n1max, n1, n2, t, w) |
| do i = 0, n2 - 1 |
| t(1, i) = 0.5d0 * (t(0, i) - t(1, i)) |
| t(0, i) = t(0, i) - t(1, i) |
| end do |
| if (n1 .gt. 4) then |
| call rftfrow(n1max, n1, n2, t, nc, w(nw)) |
| call bitrv2row(n1max, n1, n2, ip(2), t) |
| end if |
| call cftfrow(n1max, n1, n2, t, w) |
| do i = 0, n2h - 1 |
| ix = 2 * i |
| ic = n2 - 1 - i |
| do j = 0, n1h - 1 |
| jx = 2 * j |
| jc = n1 - 1 - j |
| a(jx, ix) = t(j, i) |
| a(jx + 1, ix) = t(jc, i) |
| a(jx, ix + 1) = t(j, ic) |
| a(jx + 1, ix + 1) = t(jc, ic) |
| end do |
| end do |
| else |
| do i = 0, n2h - 1 |
| ix = 2 * i |
| ic = n2 - 1 - i |
| do j = 0, n1h - 1 |
| jx = 2 * j |
| jc = n1 - 1 - j |
| t(j, i) = a(jx, ix) |
| t(jc, i) = a(jx + 1, ix) |
| t(j, ic) = a(jx, ix + 1) |
| t(jc, ic) = a(jx + 1, ix + 1) |
| end do |
| end do |
| if (n1 .gt. 4) then |
| call bitrv2row(n1max, n1, n2, ip(2), t) |
| end if |
| call cftbrow(n1max, n1, n2, t, w) |
| if (n1 .gt. 4) then |
| call rftbrow(n1max, n1, n2, t, nc, w(nw)) |
| end if |
| do i = 0, n2 - 1 |
| xi = t(0, i) - t(1, i) |
| t(0, i) = t(0, i) + t(1, i) |
| t(1, i) = xi |
| end do |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), t) |
| end if |
| call cftbcol(n1max, n1, n2, t, w) |
| call dctbsub(n1max, n1, n2, t, nc, w(nw)) |
| do i = 0, n2 - 1 |
| do j = 1, n1h - 1 |
| jx = 2 * j |
| a(j, i) = t(jx, i) |
| a(n1 - j, i) = t(jx + 1, i) |
| end do |
| end do |
| a(0, 0) = t(0, 0) |
| a(n1h, 0) = t(1, 0) |
| a(0, n2h) = t(0, n2h) |
| a(n1h, n2h) = t(1, n2h) |
| do i = 1, n2h - 1 |
| ic = n2 - i |
| a(0, i) = t(0, i) |
| a(n1h, i) = t(1, ic) |
| a(0, ic) = t(1, i) |
| a(n1h, ic) = t(0, ic) |
| end do |
| end if |
| end |
| ! |
| subroutine ddst2d(n1max, n1, n2, isgn, a, t, ip, w) |
| integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, n1h, n2h, |
| & i, ix, ic, j, jx, jc |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), |
| & t(0 : n1max - 1, 0 : n2 - 1), w(0 : *), xi |
| n = max(n1, 2 * n2) |
| nw = ip(0) |
| if (n .gt. 4 * nw) then |
| nw = n / 4 |
| call makewt(nw, ip, w) |
| end if |
| nc = ip(1) |
| if (n1 .gt. nc .or. n2 .gt. nc) then |
| nc = max(n1, n2) |
| call makect(nc, ip, w(nw)) |
| end if |
| n1h = n1 / 2 |
| n2h = n2 / 2 |
| if (isgn .ge. 0) then |
| do i = 0, n2 - 1 |
| do j = 1, n1h - 1 |
| jx = 2 * j |
| t(jx, i) = a(j, i) |
| t(jx + 1, i) = a(n1 - j, i) |
| end do |
| end do |
| t(0, 0) = a(0, 0) |
| t(1, 0) = a(n1h, 0) |
| t(0, n2h) = a(0, n2h) |
| t(1, n2h) = a(n1h, n2h) |
| do i = 1, n2h - 1 |
| ic = n2 - i |
| t(0, i) = a(0, i) |
| t(1, ic) = a(n1h, i) |
| t(1, i) = a(0, ic) |
| t(0, ic) = a(n1h, ic) |
| end do |
| call dstfsub(n1max, n1, n2, t, nc, w(nw)) |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), t) |
| end if |
| call cftfcol(n1max, n1, n2, t, w) |
| do i = 0, n2 - 1 |
| t(1, i) = 0.5d0 * (t(0, i) - t(1, i)) |
| t(0, i) = t(0, i) - t(1, i) |
| end do |
| if (n1 .gt. 4) then |
| call rftfrow(n1max, n1, n2, t, nc, w(nw)) |
| call bitrv2row(n1max, n1, n2, ip(2), t) |
| end if |
| call cftfrow(n1max, n1, n2, t, w) |
| do i = 0, n2h - 1 |
| ix = 2 * i |
| ic = n2 - 1 - i |
| do j = 0, n1h - 1 |
| jx = 2 * j |
| jc = n1 - 1 - j |
| a(jx, ix) = t(j, i) |
| a(jx + 1, ix) = -t(jc, i) |
| a(jx, ix + 1) = -t(j, ic) |
| a(jx + 1, ix + 1) = t(jc, ic) |
| end do |
| end do |
| else |
| do i = 0, n2h - 1 |
| ix = 2 * i |
| ic = n2 - 1 - i |
| do j = 0, n1h - 1 |
| jx = 2 * j |
| jc = n1 - 1 - j |
| t(j, i) = a(jx, ix) |
| t(jc, i) = -a(jx + 1, ix) |
| t(j, ic) = -a(jx, ix + 1) |
| t(jc, ic) = a(jx + 1, ix + 1) |
| end do |
| end do |
| if (n1 .gt. 4) then |
| call bitrv2row(n1max, n1, n2, ip(2), t) |
| end if |
| call cftbrow(n1max, n1, n2, t, w) |
| if (n1 .gt. 4) then |
| call rftbrow(n1max, n1, n2, t, nc, w(nw)) |
| end if |
| do i = 0, n2 - 1 |
| xi = t(0, i) - t(1, i) |
| t(0, i) = t(0, i) + t(1, i) |
| t(1, i) = xi |
| end do |
| if (n2 .gt. 2) then |
| call bitrv2col(n1max, n1, n2, ip(2), t) |
| end if |
| call cftbcol(n1max, n1, n2, t, w) |
| call dstbsub(n1max, n1, n2, t, nc, w(nw)) |
| do i = 0, n2 - 1 |
| do j = 1, n1h - 1 |
| jx = 2 * j |
| a(j, i) = t(jx, i) |
| a(n1 - j, i) = t(jx + 1, i) |
| end do |
| end do |
| a(0, 0) = t(0, 0) |
| a(n1h, 0) = t(1, 0) |
| a(0, n2h) = t(0, n2h) |
| a(n1h, n2h) = t(1, n2h) |
| do i = 1, n2h - 1 |
| ic = n2 - i |
| a(0, i) = t(0, i) |
| a(n1h, i) = t(1, ic) |
| a(0, ic) = t(1, i) |
| a(n1h, ic) = t(0, ic) |
| end do |
| end if |
| end |
| ! |
| ! -------- initializing routines -------- |
| ! |
| subroutine makewt(nw, ip, w) |
| integer nw, ip(0 : *), nwh, j |
| real*8 w(0 : nw - 1), delta, x, y |
| ip(0) = nw |
| ip(1) = 1 |
| if (nw .gt. 2) then |
| nwh = nw / 2 |
| delta = atan(1.0d0) / nwh |
| w(0) = 1 |
| w(1) = 0 |
| w(nwh) = cos(delta * nwh) |
| w(nwh + 1) = w(nwh) |
| do j = 2, nwh - 2, 2 |
| x = cos(delta * j) |
| y = sin(delta * j) |
| w(j) = x |
| w(j + 1) = y |
| w(nw - j) = y |
| w(nw - j + 1) = x |
| end do |
| call bitrv2(nw, ip(2), w) |
| end if |
| end |
| ! |
| subroutine makect(nc, ip, c) |
| integer nc, ip(0 : *), nch, j |
| real*8 c(0 : nc - 1), delta |
| ip(1) = nc |
| if (nc .gt. 1) then |
| nch = nc / 2 |
| delta = atan(1.0d0) / nch |
| c(0) = 0.5d0 |
| c(nch) = 0.5d0 * cos(delta * nch) |
| do j = 1, nch - 1 |
| c(j) = 0.5d0 * cos(delta * j) |
| c(nc - j) = 0.5d0 * sin(delta * j) |
| end do |
| end if |
| end |
| ! |
| ! -------- child routines -------- |
| ! |
| subroutine bitrv2(n, ip, a) |
| integer n, ip(0 : *), j, j1, k, k1, l, m, m2 |
| real*8 a(0 : n - 1), xr, xi |
| ip(0) = 0 |
| l = n |
| m = 1 |
| do while (4 * m .lt. l) |
| l = l / 2 |
| do j = 0, m - 1 |
| ip(m + j) = ip(j) + l |
| end do |
| m = m * 2 |
| end do |
| if (4 * m .gt. l) then |
| do k = 1, m - 1 |
| do j = 0, k - 1 |
| j1 = 2 * j + ip(k) |
| k1 = 2 * k + ip(j) |
| xr = a(j1) |
| xi = a(j1 + 1) |
| a(j1) = a(k1) |
| a(j1 + 1) = a(k1 + 1) |
| a(k1) = xr |
| a(k1 + 1) = xi |
| end do |
| end do |
| else |
| m2 = 2 * m |
| do k = 1, m - 1 |
| do j = 0, k - 1 |
| j1 = 2 * j + ip(k) |
| k1 = 2 * k + ip(j) |
| xr = a(j1) |
| xi = a(j1 + 1) |
| a(j1) = a(k1) |
| a(j1 + 1) = a(k1 + 1) |
| a(k1) = xr |
| a(k1 + 1) = xi |
| j1 = j1 + m2 |
| k1 = k1 + m2 |
| xr = a(j1) |
| xi = a(j1 + 1) |
| a(j1) = a(k1) |
| a(j1 + 1) = a(k1 + 1) |
| a(k1) = xr |
| a(k1 + 1) = xi |
| end do |
| end do |
| end if |
| end |
| ! |
| subroutine bitrv2row(n1max, n, n2, ip, a) |
| integer n1max, n, n2, ip(0 : *), i, j, j1, k, k1, l, m, m2 |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), xr, xi |
| ip(0) = 0 |
| l = n |
| m = 1 |
| do while (4 * m .lt. l) |
| l = l / 2 |
| do j = 0, m - 1 |
| ip(m + j) = ip(j) + l |
| end do |
| m = m * 2 |
| end do |
| if (4 * m .gt. l) then |
| do i = 0, n2 - 1 |
| do k = 1, m - 1 |
| do j = 0, k - 1 |
| j1 = 2 * j + ip(k) |
| k1 = 2 * k + ip(j) |
| xr = a(j1, i) |
| xi = a(j1 + 1, i) |
| a(j1, i) = a(k1, i) |
| a(j1 + 1, i) = a(k1 + 1, i) |
| a(k1, i) = xr |
| a(k1 + 1, i) = xi |
| end do |
| end do |
| end do |
| else |
| m2 = 2 * m |
| do i = 0, n2 - 1 |
| do k = 1, m - 1 |
| do j = 0, k - 1 |
| j1 = 2 * j + ip(k) |
| k1 = 2 * k + ip(j) |
| xr = a(j1, i) |
| xi = a(j1 + 1, i) |
| a(j1, i) = a(k1, i) |
| a(j1 + 1, i) = a(k1 + 1, i) |
| a(k1, i) = xr |
| a(k1 + 1, i) = xi |
| j1 = j1 + m2 |
| k1 = k1 + m2 |
| xr = a(j1, i) |
| xi = a(j1 + 1, i) |
| a(j1, i) = a(k1, i) |
| a(j1 + 1, i) = a(k1 + 1, i) |
| a(k1, i) = xr |
| a(k1 + 1, i) = xi |
| end do |
| end do |
| end do |
| end if |
| end |
| ! |
| subroutine bitrv2col(n1max, n1, n, ip, a) |
| integer n1max, n1, n, ip(0 : *), i, j, j1, k, k1, l, m |
| real*8 a(0 : n1max - 1, 0 : n - 1), xr, xi |
| ip(0) = 0 |
| l = n |
| m = 1 |
| do while (2 * m .lt. l) |
| l = l / 2 |
| do j = 0, m - 1 |
| ip(m + j) = ip(j) + l |
| end do |
| m = m * 2 |
| end do |
| if (2 * m .gt. l) then |
| do k = 1, m - 1 |
| do j = 0, k - 1 |
| j1 = j + ip(k) |
| k1 = k + ip(j) |
| do i = 0, n1 - 2, 2 |
| xr = a(i, j1) |
| xi = a(i + 1, j1) |
| a(i, j1) = a(i, k1) |
| a(i + 1, j1) = a(i + 1, k1) |
| a(i, k1) = xr |
| a(i + 1, k1) = xi |
| end do |
| end do |
| end do |
| else |
| do k = 1, m - 1 |
| do j = 0, k - 1 |
| j1 = j + ip(k) |
| k1 = k + ip(j) |
| do i = 0, n1 - 2, 2 |
| xr = a(i, j1) |
| xi = a(i + 1, j1) |
| a(i, j1) = a(i, k1) |
| a(i + 1, j1) = a(i + 1, k1) |
| a(i, k1) = xr |
| a(i + 1, k1) = xi |
| end do |
| j1 = j1 + m |
| k1 = k1 + m |
| do i = 0, n1 - 2, 2 |
| xr = a(i, j1) |
| xi = a(i + 1, j1) |
| a(i, j1) = a(i, k1) |
| a(i + 1, j1) = a(i + 1, k1) |
| a(i, k1) = xr |
| a(i + 1, k1) = xi |
| end do |
| end do |
| end do |
| end if |
| end |
| ! |
| subroutine cftbrow(n1max, n, n2, a, w) |
| integer n1max, n, n2, i, j, j1, j2, j3, k, k1, ks, l, m |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *) |
| real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i |
| real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i |
| do i = 0, n2 - 1 |
| l = 2 |
| do while (2 * l .lt. n) |
| m = 4 * l |
| do j = 0, l - 2, 2 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| x0r = a(j, i) + a(j1, i) |
| x0i = a(j + 1, i) + a(j1 + 1, i) |
| x1r = a(j, i) - a(j1, i) |
| x1i = a(j + 1, i) - a(j1 + 1, i) |
| x2r = a(j2, i) + a(j3, i) |
| x2i = a(j2 + 1, i) + a(j3 + 1, i) |
| x3r = a(j2, i) - a(j3, i) |
| x3i = a(j2 + 1, i) - a(j3 + 1, i) |
| a(j, i) = x0r + x2r |
| a(j + 1, i) = x0i + x2i |
| a(j2, i) = x0r - x2r |
| a(j2 + 1, i) = x0i - x2i |
| a(j1, i) = x1r - x3i |
| a(j1 + 1, i) = x1i + x3r |
| a(j3, i) = x1r + x3i |
| a(j3 + 1, i) = x1i - x3r |
| end do |
| if (m .lt. n) then |
| wk1r = w(2) |
| do j = m, l + m - 2, 2 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| x0r = a(j, i) + a(j1, i) |
| x0i = a(j + 1, i) + a(j1 + 1, i) |
| x1r = a(j, i) - a(j1, i) |
| x1i = a(j + 1, i) - a(j1 + 1, i) |
| x2r = a(j2, i) + a(j3, i) |
| x2i = a(j2 + 1, i) + a(j3 + 1, i) |
| x3r = a(j2, i) - a(j3, i) |
| x3i = a(j2 + 1, i) - a(j3 + 1, i) |
| a(j, i) = x0r + x2r |
| a(j + 1, i) = x0i + x2i |
| a(j2, i) = x2i - x0i |
| a(j2 + 1, i) = x0r - x2r |
| x0r = x1r - x3i |
| x0i = x1i + x3r |
| a(j1, i) = wk1r * (x0r - x0i) |
| a(j1 + 1, i) = wk1r * (x0r + x0i) |
| x0r = x3i + x1r |
| x0i = x3r - x1i |
| a(j3, i) = wk1r * (x0i - x0r) |
| a(j3 + 1, i) = wk1r * (x0i + x0r) |
| end do |
| k1 = 1 |
| ks = -1 |
| do k = 2 * m, n - m, m |
| k1 = k1 + 1 |
| ks = -ks |
| wk1r = w(2 * k1) |
| wk1i = w(2 * k1 + 1) |
| wk2r = ks * w(k1) |
| wk2i = w(k1 + ks) |
| wk3r = wk1r - 2 * wk2i * wk1i |
| wk3i = 2 * wk2i * wk1r - wk1i |
| do j = k, l + k - 2, 2 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| x0r = a(j, i) + a(j1, i) |
| x0i = a(j + 1, i) + a(j1 + 1, i) |
| x1r = a(j, i) - a(j1, i) |
| x1i = a(j + 1, i) - a(j1 + 1, i) |
| x2r = a(j2, i) + a(j3, i) |
| x2i = a(j2 + 1, i) + a(j3 + 1, i) |
| x3r = a(j2, i) - a(j3, i) |
| x3i = a(j2 + 1, i) - a(j3 + 1, i) |
| a(j, i) = x0r + x2r |
| a(j + 1, i) = x0i + x2i |
| x0r = x0r - x2r |
| x0i = x0i - x2i |
| a(j2, i) = wk2r * x0r - wk2i * x0i |
| a(j2 + 1, i) = wk2r * x0i + wk2i * x0r |
| x0r = x1r - x3i |
| x0i = x1i + x3r |
| a(j1, i) = wk1r * x0r - wk1i * x0i |
| a(j1 + 1, i) = wk1r * x0i + wk1i * x0r |
| x0r = x1r + x3i |
| x0i = x1i - x3r |
| a(j3, i) = wk3r * x0r - wk3i * x0i |
| a(j3 + 1, i) = wk3r * x0i + wk3i * x0r |
| end do |
| end do |
| end if |
| l = m |
| end do |
| if (l .lt. n) then |
| do j = 0, l - 2, 2 |
| j1 = j + l |
| x0r = a(j, i) - a(j1, i) |
| x0i = a(j + 1, i) - a(j1 + 1, i) |
| a(j, i) = a(j, i) + a(j1, i) |
| a(j + 1, i) = a(j + 1, i) + a(j1 + 1, i) |
| a(j1, i) = x0r |
| a(j1 + 1, i) = x0i |
| end do |
| end if |
| end do |
| end |
| ! |
| subroutine cftbcol(n1max, n1, n, a, w) |
| integer n1max, n1, n, i, j, j1, j2, j3, k, k1, ks, l, m |
| real*8 a(0 : n1max - 1, 0 : n - 1), w(0 : *) |
| real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i |
| real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i |
| l = 1 |
| do while (2 * l .lt. n) |
| m = 4 * l |
| do j = 0, l - 1 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) + a(i, j1) |
| x0i = a(i + 1, j) + a(i + 1, j1) |
| x1r = a(i, j) - a(i, j1) |
| x1i = a(i + 1, j) - a(i + 1, j1) |
| x2r = a(i, j2) + a(i, j3) |
| x2i = a(i + 1, j2) + a(i + 1, j3) |
| x3r = a(i, j2) - a(i, j3) |
| x3i = a(i + 1, j2) - a(i + 1, j3) |
| a(i, j) = x0r + x2r |
| a(i + 1, j) = x0i + x2i |
| a(i, j2) = x0r - x2r |
| a(i + 1, j2) = x0i - x2i |
| a(i, j1) = x1r - x3i |
| a(i + 1, j1) = x1i + x3r |
| a(i, j3) = x1r + x3i |
| a(i + 1, j3) = x1i - x3r |
| end do |
| end do |
| if (m .lt. n) then |
| wk1r = w(2) |
| do j = m, l + m - 1 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) + a(i, j1) |
| x0i = a(i + 1, j) + a(i + 1, j1) |
| x1r = a(i, j) - a(i, j1) |
| x1i = a(i + 1, j) - a(i + 1, j1) |
| x2r = a(i, j2) + a(i, j3) |
| x2i = a(i + 1, j2) + a(i + 1, j3) |
| x3r = a(i, j2) - a(i, j3) |
| x3i = a(i + 1, j2) - a(i + 1, j3) |
| a(i, j) = x0r + x2r |
| a(i + 1, j) = x0i + x2i |
| a(i, j2) = x2i - x0i |
| a(i + 1, j2) = x0r - x2r |
| x0r = x1r - x3i |
| x0i = x1i + x3r |
| a(i, j1) = wk1r * (x0r - x0i) |
| a(i + 1, j1) = wk1r * (x0r + x0i) |
| x0r = x3i + x1r |
| x0i = x3r - x1i |
| a(i, j3) = wk1r * (x0i - x0r) |
| a(i + 1, j3) = wk1r * (x0i + x0r) |
| end do |
| end do |
| k1 = 1 |
| ks = -1 |
| do k = 2 * m, n - m, m |
| k1 = k1 + 1 |
| ks = -ks |
| wk1r = w(2 * k1) |
| wk1i = w(2 * k1 + 1) |
| wk2r = ks * w(k1) |
| wk2i = w(k1 + ks) |
| wk3r = wk1r - 2 * wk2i * wk1i |
| wk3i = 2 * wk2i * wk1r - wk1i |
| do j = k, l + k - 1 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) + a(i, j1) |
| x0i = a(i + 1, j) + a(i + 1, j1) |
| x1r = a(i, j) - a(i, j1) |
| x1i = a(i + 1, j) - a(i + 1, j1) |
| x2r = a(i, j2) + a(i, j3) |
| x2i = a(i + 1, j2) + a(i + 1, j3) |
| x3r = a(i, j2) - a(i, j3) |
| x3i = a(i + 1, j2) - a(i + 1, j3) |
| a(i, j) = x0r + x2r |
| a(i + 1, j) = x0i + x2i |
| x0r = x0r - x2r |
| x0i = x0i - x2i |
| a(i, j2) = wk2r * x0r - wk2i * x0i |
| a(i + 1, j2) = wk2r * x0i + wk2i * x0r |
| x0r = x1r - x3i |
| x0i = x1i + x3r |
| a(i, j1) = wk1r * x0r - wk1i * x0i |
| a(i + 1, j1) = wk1r * x0i + wk1i * x0r |
| x0r = x1r + x3i |
| x0i = x1i - x3r |
| a(i, j3) = wk3r * x0r - wk3i * x0i |
| a(i + 1, j3) = wk3r * x0i + wk3i * x0r |
| end do |
| end do |
| end do |
| end if |
| l = m |
| end do |
| if (l .lt. n) then |
| do j = 0, l - 1 |
| j1 = j + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) - a(i, j1) |
| x0i = a(i + 1, j) - a(i + 1, j1) |
| a(i, j) = a(i, j) + a(i, j1) |
| a(i + 1, j) = a(i + 1, j) + a(i + 1, j1) |
| a(i, j1) = x0r |
| a(i + 1, j1) = x0i |
| end do |
| end do |
| end if |
| end |
| ! |
| subroutine cftfrow(n1max, n, n2, a, w) |
| integer n1max, n, n2, i, j, j1, j2, j3, k, k1, ks, l, m |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *) |
| real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i |
| real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i |
| do i = 0, n2 - 1 |
| l = 2 |
| do while (2 * l .lt. n) |
| m = 4 * l |
| do j = 0, l - 2, 2 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| x0r = a(j, i) + a(j1, i) |
| x0i = a(j + 1, i) + a(j1 + 1, i) |
| x1r = a(j, i) - a(j1, i) |
| x1i = a(j + 1, i) - a(j1 + 1, i) |
| x2r = a(j2, i) + a(j3, i) |
| x2i = a(j2 + 1, i) + a(j3 + 1, i) |
| x3r = a(j2, i) - a(j3, i) |
| x3i = a(j2 + 1, i) - a(j3 + 1, i) |
| a(j, i) = x0r + x2r |
| a(j + 1, i) = x0i + x2i |
| a(j2, i) = x0r - x2r |
| a(j2 + 1, i) = x0i - x2i |
| a(j1, i) = x1r + x3i |
| a(j1 + 1, i) = x1i - x3r |
| a(j3, i) = x1r - x3i |
| a(j3 + 1, i) = x1i + x3r |
| end do |
| if (m .lt. n) then |
| wk1r = w(2) |
| do j = m, l + m - 2, 2 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| x0r = a(j, i) + a(j1, i) |
| x0i = a(j + 1, i) + a(j1 + 1, i) |
| x1r = a(j, i) - a(j1, i) |
| x1i = a(j + 1, i) - a(j1 + 1, i) |
| x2r = a(j2, i) + a(j3, i) |
| x2i = a(j2 + 1, i) + a(j3 + 1, i) |
| x3r = a(j2, i) - a(j3, i) |
| x3i = a(j2 + 1, i) - a(j3 + 1, i) |
| a(j, i) = x0r + x2r |
| a(j + 1, i) = x0i + x2i |
| a(j2, i) = x0i - x2i |
| a(j2 + 1, i) = x2r - x0r |
| x0r = x1r + x3i |
| x0i = x1i - x3r |
| a(j1, i) = wk1r * (x0i + x0r) |
| a(j1 + 1, i) = wk1r * (x0i - x0r) |
| x0r = x3i - x1r |
| x0i = x3r + x1i |
| a(j3, i) = wk1r * (x0r + x0i) |
| a(j3 + 1, i) = wk1r * (x0r - x0i) |
| end do |
| k1 = 1 |
| ks = -1 |
| do k = 2 * m, n - m, m |
| k1 = k1 + 1 |
| ks = -ks |
| wk1r = w(2 * k1) |
| wk1i = w(2 * k1 + 1) |
| wk2r = ks * w(k1) |
| wk2i = w(k1 + ks) |
| wk3r = wk1r - 2 * wk2i * wk1i |
| wk3i = 2 * wk2i * wk1r - wk1i |
| do j = k, l + k - 2, 2 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| x0r = a(j, i) + a(j1, i) |
| x0i = a(j + 1, i) + a(j1 + 1, i) |
| x1r = a(j, i) - a(j1, i) |
| x1i = a(j + 1, i) - a(j1 + 1, i) |
| x2r = a(j2, i) + a(j3, i) |
| x2i = a(j2 + 1, i) + a(j3 + 1, i) |
| x3r = a(j2, i) - a(j3, i) |
| x3i = a(j2 + 1, i) - a(j3 + 1, i) |
| a(j, i) = x0r + x2r |
| a(j + 1, i) = x0i + x2i |
| x0r = x0r - x2r |
| x0i = x0i - x2i |
| a(j2, i) = wk2r * x0r + wk2i * x0i |
| a(j2 + 1, i) = wk2r * x0i - wk2i * x0r |
| x0r = x1r + x3i |
| x0i = x1i - x3r |
| a(j1, i) = wk1r * x0r + wk1i * x0i |
| a(j1 + 1, i) = wk1r * x0i - wk1i * x0r |
| x0r = x1r - x3i |
| x0i = x1i + x3r |
| a(j3, i) = wk3r * x0r + wk3i * x0i |
| a(j3 + 1, i) = wk3r * x0i - wk3i * x0r |
| end do |
| end do |
| end if |
| l = m |
| end do |
| if (l .lt. n) then |
| do j = 0, l - 2, 2 |
| j1 = j + l |
| x0r = a(j, i) - a(j1, i) |
| x0i = a(j + 1, i) - a(j1 + 1, i) |
| a(j, i) = a(j, i) + a(j1, i) |
| a(j + 1, i) = a(j + 1, i) + a(j1 + 1, i) |
| a(j1, i) = x0r |
| a(j1 + 1, i) = x0i |
| end do |
| end if |
| end do |
| end |
| ! |
| subroutine cftfcol(n1max, n1, n, a, w) |
| integer n1max, n1, n, i, j, j1, j2, j3, k, k1, ks, l, m |
| real*8 a(0 : n1max - 1, 0 : n - 1), w(0 : *) |
| real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i |
| real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i |
| l = 1 |
| do while (2 * l .lt. n) |
| m = 4 * l |
| do j = 0, l - 1 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) + a(i, j1) |
| x0i = a(i + 1, j) + a(i + 1, j1) |
| x1r = a(i, j) - a(i, j1) |
| x1i = a(i + 1, j) - a(i + 1, j1) |
| x2r = a(i, j2) + a(i, j3) |
| x2i = a(i + 1, j2) + a(i + 1, j3) |
| x3r = a(i, j2) - a(i, j3) |
| x3i = a(i + 1, j2) - a(i + 1, j3) |
| a(i, j) = x0r + x2r |
| a(i + 1, j) = x0i + x2i |
| a(i, j2) = x0r - x2r |
| a(i + 1, j2) = x0i - x2i |
| a(i, j1) = x1r + x3i |
| a(i + 1, j1) = x1i - x3r |
| a(i, j3) = x1r - x3i |
| a(i + 1, j3) = x1i + x3r |
| end do |
| end do |
| if (m .lt. n) then |
| wk1r = w(2) |
| do j = m, l + m - 1 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) + a(i, j1) |
| x0i = a(i + 1, j) + a(i + 1, j1) |
| x1r = a(i, j) - a(i, j1) |
| x1i = a(i + 1, j) - a(i + 1, j1) |
| x2r = a(i, j2) + a(i, j3) |
| x2i = a(i + 1, j2) + a(i + 1, j3) |
| x3r = a(i, j2) - a(i, j3) |
| x3i = a(i + 1, j2) - a(i + 1, j3) |
| a(i, j) = x0r + x2r |
| a(i + 1, j) = x0i + x2i |
| a(i, j2) = x0i - x2i |
| a(i + 1, j2) = x2r - x0r |
| x0r = x1r + x3i |
| x0i = x1i - x3r |
| a(i, j1) = wk1r * (x0i + x0r) |
| a(i + 1, j1) = wk1r * (x0i - x0r) |
| x0r = x3i - x1r |
| x0i = x3r + x1i |
| a(i, j3) = wk1r * (x0r + x0i) |
| a(i + 1, j3) = wk1r * (x0r - x0i) |
| end do |
| end do |
| k1 = 1 |
| ks = -1 |
| do k = 2 * m, n - m, m |
| k1 = k1 + 1 |
| ks = -ks |
| wk1r = w(2 * k1) |
| wk1i = w(2 * k1 + 1) |
| wk2r = ks * w(k1) |
| wk2i = w(k1 + ks) |
| wk3r = wk1r - 2 * wk2i * wk1i |
| wk3i = 2 * wk2i * wk1r - wk1i |
| do j = k, l + k - 1 |
| j1 = j + l |
| j2 = j1 + l |
| j3 = j2 + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) + a(i, j1) |
| x0i = a(i + 1, j) + a(i + 1, j1) |
| x1r = a(i, j) - a(i, j1) |
| x1i = a(i + 1, j) - a(i + 1, j1) |
| x2r = a(i, j2) + a(i, j3) |
| x2i = a(i + 1, j2) + a(i + 1, j3) |
| x3r = a(i, j2) - a(i, j3) |
| x3i = a(i + 1, j2) - a(i + 1, j3) |
| a(i, j) = x0r + x2r |
| a(i + 1, j) = x0i + x2i |
| x0r = x0r - x2r |
| x0i = x0i - x2i |
| a(i, j2) = wk2r * x0r + wk2i * x0i |
| a(i + 1, j2) = wk2r * x0i - wk2i * x0r |
| x0r = x1r + x3i |
| x0i = x1i - x3r |
| a(i, j1) = wk1r * x0r + wk1i * x0i |
| a(i + 1, j1) = wk1r * x0i - wk1i * x0r |
| x0r = x1r - x3i |
| x0i = x1i + x3r |
| a(i, j3) = wk3r * x0r + wk3i * x0i |
| a(i + 1, j3) = wk3r * x0i - wk3i * x0r |
| end do |
| end do |
| end do |
| end if |
| l = m |
| end do |
| if (l .lt. n) then |
| do j = 0, l - 1 |
| j1 = j + l |
| do i = 0, n1 - 2, 2 |
| x0r = a(i, j) - a(i, j1) |
| x0i = a(i + 1, j) - a(i + 1, j1) |
| a(i, j) = a(i, j) + a(i, j1) |
| a(i + 1, j) = a(i + 1, j) + a(i + 1, j1) |
| a(i, j1) = x0r |
| a(i + 1, j1) = x0i |
| end do |
| end do |
| end if |
| end |
| ! |
| subroutine rftbrow(n1max, n, n2, a, nc, c) |
| integer n1max, n, n2, nc, i, j, k, kk, ks |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), |
| & wkr, wki, xr, xi, yr, yi |
| ks = 4 * nc / n |
| do i = 0, n2 - 1 |
| kk = 0 |
| do k = n / 2 - 2, 2, -2 |
| j = n - k |
| kk = kk + ks |
| wkr = 0.5d0 - c(kk) |
| wki = c(nc - kk) |
| xr = a(k, i) - a(j, i) |
| xi = a(k + 1, i) + a(j + 1, i) |
| yr = wkr * xr - wki * xi |
| yi = wkr * xi + wki * xr |
| a(k, i) = a(k, i) - yr |
| a(k + 1, i) = a(k + 1, i) - yi |
| a(j, i) = a(j, i) + yr |
| a(j + 1, i) = a(j + 1, i) - yi |
| end do |
| end do |
| end |
| ! |
| subroutine rftfrow(n1max, n, n2, a, nc, c) |
| integer n1max, n, n2, nc, i, j, k, kk, ks |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), |
| & wkr, wki, xr, xi, yr, yi |
| ks = 4 * nc / n |
| do i = 0, n2 - 1 |
| kk = 0 |
| do k = n / 2 - 2, 2, -2 |
| j = n - k |
| kk = kk + ks |
| wkr = 0.5d0 - c(kk) |
| wki = c(nc - kk) |
| xr = a(k, i) - a(j, i) |
| xi = a(k + 1, i) + a(j + 1, i) |
| yr = wkr * xr + wki * xi |
| yi = wkr * xi - wki * xr |
| a(k, i) = a(k, i) - yr |
| a(k + 1, i) = a(k + 1, i) - yi |
| a(j, i) = a(j, i) + yr |
| a(j + 1, i) = a(j + 1, i) - yi |
| end do |
| end do |
| end |
| ! |
| subroutine dctbsub(n1max, n1, n2, a, nc, c) |
| integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, |
| & k1, k2 |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), |
| & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i |
| ks1 = nc / n1 |
| ks2 = nc / n2 |
| n2h = n2 / 2 |
| kk2 = ks2 |
| do k2 = 1, n2h - 1 |
| j2 = n2 - k2 |
| w2r = 2 * c(kk2) |
| w2i = 2 * c(nc - kk2) |
| kk2 = kk2 + ks2 |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| x0r = w2r * c(kk1) |
| x0i = w2i * c(kk1) |
| x1r = w2r * c(nc - kk1) |
| x1i = w2i * c(nc - kk1) |
| wkr = x0r - x1i |
| wki = x0i + x1r |
| wji = x0r + x1i |
| wjr = x0i - x1r |
| kk1 = kk1 + ks1 |
| x0r = wkr * a(k1, k2) - wki * a(k1 + 1, k2) |
| x0i = wkr * a(k1 + 1, k2) + wki * a(k1, k2) |
| x1r = wjr * a(k1, j2) - wji * a(k1 + 1, j2) |
| x1i = wjr * a(k1 + 1, j2) + wji * a(k1, j2) |
| a(k1, k2) = x0r + x1i |
| a(k1 + 1, k2) = x0i - x1r |
| a(k1, j2) = x1r + x0i |
| a(k1 + 1, j2) = x1i - x0r |
| end do |
| wkr = w2r * 0.5d0 |
| wki = w2i * 0.5d0 |
| wjr = w2r * c(kk1) |
| wji = w2i * c(kk1) |
| x0r = a(0, k2) + a(0, j2) |
| x0i = a(1, k2) - a(1, j2) |
| x1r = a(0, k2) - a(0, j2) |
| x1i = a(1, k2) + a(1, j2) |
| a(0, k2) = wkr * x0r - wki * x0i |
| a(1, k2) = wkr * x0i + wki * x0r |
| a(0, j2) = -wjr * x1r + wji * x1i |
| a(1, j2) = wjr * x1i + wji * x1r |
| end do |
| w2r = 2 * c(kk2) |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| wkr = 2 * c(kk1) |
| wki = 2 * c(nc - kk1) |
| wjr = w2r * wkr |
| wji = w2r * wki |
| kk1 = kk1 + ks1 |
| x0i = wkr * a(k1 + 1, 0) + wki * a(k1, 0) |
| a(k1, 0) = wkr * a(k1, 0) - wki * a(k1 + 1, 0) |
| a(k1 + 1, 0) = x0i |
| x0i = wjr * a(k1 + 1, n2h) + wji * a(k1, n2h) |
| a(k1, n2h) = wjr * a(k1, n2h) - wji * a(k1 + 1, n2h) |
| a(k1 + 1, n2h) = x0i |
| end do |
| a(1, 0) = a(1, 0) * w2r |
| a(0, n2h) = a(0, n2h) * w2r |
| a(1, n2h) = a(1, n2h) * 0.5d0 |
| end |
| ! |
| subroutine dctfsub(n1max, n1, n2, a, nc, c) |
| integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, |
| & k1, k2 |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), |
| & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i |
| ks1 = nc / n1 |
| ks2 = nc / n2 |
| n2h = n2 / 2 |
| kk2 = ks2 |
| do k2 = 1, n2h - 1 |
| j2 = n2 - k2 |
| w2r = 2 * c(kk2) |
| w2i = 2 * c(nc - kk2) |
| kk2 = kk2 + ks2 |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| x0r = w2r * c(kk1) |
| x0i = w2i * c(kk1) |
| x1r = w2r * c(nc - kk1) |
| x1i = w2i * c(nc - kk1) |
| wkr = x0r - x1i |
| wki = x0i + x1r |
| wji = x0r + x1i |
| wjr = x0i - x1r |
| kk1 = kk1 + ks1 |
| x0r = a(k1, k2) - a(k1 + 1, j2) |
| x0i = a(k1, j2) + a(k1 + 1, k2) |
| x1r = a(k1, j2) - a(k1 + 1, k2) |
| x1i = a(k1, k2) + a(k1 + 1, j2) |
| a(k1, k2) = wkr * x0r + wki * x0i |
| a(k1 + 1, k2) = wkr * x0i - wki * x0r |
| a(k1, j2) = wjr * x1r + wji * x1i |
| a(k1 + 1, j2) = wjr * x1i - wji * x1r |
| end do |
| x0r = 2 * c(kk1) |
| wjr = x0r * w2r |
| wji = x0r * w2i |
| x0r = w2r * a(0, k2) + w2i * a(1, k2) |
| x0i = w2r * a(1, k2) - w2i * a(0, k2) |
| x1r = -wjr * a(0, j2) + wji * a(1, j2) |
| x1i = wjr * a(1, j2) + wji * a(0, j2) |
| a(0, k2) = x0r + x1r |
| a(1, k2) = x1i + x0i |
| a(0, j2) = x0r - x1r |
| a(1, j2) = x1i - x0i |
| end do |
| w2r = 2 * c(kk2) |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| wkr = 2 * c(kk1) |
| wki = 2 * c(nc - kk1) |
| wjr = w2r * wkr |
| wji = w2r * wki |
| kk1 = kk1 + ks1 |
| x0i = wkr * a(k1 + 1, 0) - wki * a(k1, 0) |
| a(k1, 0) = wkr * a(k1, 0) + wki * a(k1 + 1, 0) |
| a(k1 + 1, 0) = x0i |
| x0i = wjr * a(k1 + 1, n2h) - wji * a(k1, n2h) |
| a(k1, n2h) = wjr * a(k1, n2h) + wji * a(k1 + 1, n2h) |
| a(k1 + 1, n2h) = x0i |
| end do |
| w2r = w2r * 2 |
| a(0, 0) = a(0, 0) * 2 |
| a(1, 0) = a(1, 0) * w2r |
| a(0, n2h) = a(0, n2h) * w2r |
| end |
| ! |
| subroutine dstbsub(n1max, n1, n2, a, nc, c) |
| integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, |
| & k1, k2 |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), |
| & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i |
| ks1 = nc / n1 |
| ks2 = nc / n2 |
| n2h = n2 / 2 |
| kk2 = ks2 |
| do k2 = 1, n2h - 1 |
| j2 = n2 - k2 |
| w2r = 2 * c(kk2) |
| w2i = 2 * c(nc - kk2) |
| kk2 = kk2 + ks2 |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| x0r = w2r * c(kk1) |
| x0i = w2i * c(kk1) |
| x1r = w2r * c(nc - kk1) |
| x1i = w2i * c(nc - kk1) |
| wkr = x0r - x1i |
| wki = x0i + x1r |
| wji = x0r + x1i |
| wjr = x0i - x1r |
| kk1 = kk1 + ks1 |
| x0r = wkr * a(k1, k2) - wki * a(k1 + 1, k2) |
| x0i = wkr * a(k1 + 1, k2) + wki * a(k1, k2) |
| x1r = wjr * a(k1, j2) - wji * a(k1 + 1, j2) |
| x1i = wjr * a(k1 + 1, j2) + wji * a(k1, j2) |
| a(k1, k2) = x1i - x0r |
| a(k1 + 1, k2) = x1r + x0i |
| a(k1, j2) = x0i - x1r |
| a(k1 + 1, j2) = x0r + x1i |
| end do |
| wkr = w2r * 0.5d0 |
| wki = w2i * 0.5d0 |
| wjr = w2r * c(kk1) |
| wji = w2i * c(kk1) |
| x0r = a(0, k2) + a(0, j2) |
| x0i = a(1, k2) - a(1, j2) |
| x1r = a(0, k2) - a(0, j2) |
| x1i = a(1, k2) + a(1, j2) |
| a(1, k2) = wkr * x0r - wki * x0i |
| a(0, k2) = wkr * x0i + wki * x0r |
| a(1, j2) = -wjr * x1r + wji * x1i |
| a(0, j2) = wjr * x1i + wji * x1r |
| end do |
| w2r = 2 * c(kk2) |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| wkr = 2 * c(kk1) |
| wki = 2 * c(nc - kk1) |
| wjr = w2r * wkr |
| wji = w2r * wki |
| kk1 = kk1 + ks1 |
| x0i = wkr * a(k1 + 1, 0) + wki * a(k1, 0) |
| a(k1 + 1, 0) = wkr * a(k1, 0) - wki * a(k1 + 1, 0) |
| a(k1, 0) = x0i |
| x0i = wjr * a(k1 + 1, n2h) + wji * a(k1, n2h) |
| a(k1 + 1, n2h) = wjr * a(k1, n2h) - wji * a(k1 + 1, n2h) |
| a(k1, n2h) = x0i |
| end do |
| a(1, 0) = a(1, 0) * w2r |
| a(0, n2h) = a(0, n2h) * w2r |
| a(1, n2h) = a(1, n2h) * 0.5d0 |
| end |
| ! |
| subroutine dstfsub(n1max, n1, n2, a, nc, c) |
| integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, |
| & k1, k2 |
| real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), |
| & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i |
| ks1 = nc / n1 |
| ks2 = nc / n2 |
| n2h = n2 / 2 |
| kk2 = ks2 |
| do k2 = 1, n2h - 1 |
| j2 = n2 - k2 |
| w2r = 2 * c(kk2) |
| w2i = 2 * c(nc - kk2) |
| kk2 = kk2 + ks2 |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| x0r = w2r * c(kk1) |
| x0i = w2i * c(kk1) |
| x1r = w2r * c(nc - kk1) |
| x1i = w2i * c(nc - kk1) |
| wkr = x0r - x1i |
| wki = x0i + x1r |
| wji = x0r + x1i |
| wjr = x0i - x1r |
| kk1 = kk1 + ks1 |
| x0r = a(k1 + 1, j2) - a(k1, k2) |
| x0i = a(k1 + 1, k2) + a(k1, j2) |
| x1r = a(k1 + 1, k2) - a(k1, j2) |
| x1i = a(k1 + 1, j2) + a(k1, k2) |
| a(k1, k2) = wkr * x0r + wki * x0i |
| a(k1 + 1, k2) = wkr * x0i - wki * x0r |
| a(k1, j2) = wjr * x1r + wji * x1i |
| a(k1 + 1, j2) = wjr * x1i - wji * x1r |
| end do |
| x0r = 2 * c(kk1) |
| wjr = x0r * w2r |
| wji = x0r * w2i |
| x0r = w2r * a(1, k2) + w2i * a(0, k2) |
| x0i = w2r * a(0, k2) - w2i * a(1, k2) |
| x1r = -wjr * a(1, j2) + wji * a(0, j2) |
| x1i = wjr * a(0, j2) + wji * a(1, j2) |
| a(0, k2) = x0r + x1r |
| a(1, k2) = x1i + x0i |
| a(0, j2) = x0r - x1r |
| a(1, j2) = x1i - x0i |
| end do |
| w2r = 2 * c(kk2) |
| kk1 = ks1 |
| do k1 = 2, n1 - 2, 2 |
| wkr = 2 * c(kk1) |
| wki = 2 * c(nc - kk1) |
| wjr = w2r * wkr |
| wji = w2r * wki |
| kk1 = kk1 + ks1 |
| x0i = wkr * a(k1, 0) - wki * a(k1 + 1, 0) |
| a(k1, 0) = wkr * a(k1 + 1, 0) + wki * a(k1, 0) |
| a(k1 + 1, 0) = x0i |
| x0i = wjr * a(k1, n2h) - wji * a(k1 + 1, n2h) |
| a(k1, n2h) = wjr * a(k1 + 1, n2h) + wji * a(k1, n2h) |
| a(k1 + 1, n2h) = x0i |
| end do |
| w2r = w2r * 2 |
| a(0, 0) = a(0, 0) * 2 |
| a(1, 0) = a(1, 0) * w2r |
| a(0, n2h) = a(0, n2h) * w2r |
| end |
| ! |
