fft4f2d.f

2维fft程序

                  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