fft4f2d.f

2维fft程序

                      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
!