fftsg.f

2维fft程序

              yr = a(k) - a(n - k)
              yi = a(k) + a(n - k)
              a(j) = xr
              a(k) = yr
              t(j) = xi - yi
              t(k) = xi + yi
          end do
          t(mh) = a(mh) + a(n - mh)
          a(mh) = a(mh) - a(n - mh)
          call dctsub(m, a, nc, w(nw))
          if (m .gt. 4) then
              call cftfsub(m, a, ip, nw, w)
              call rftfsub(m, a, nc, w(nw))
          else if (m .eq. 4) then
              call cftfsub(m, a, ip, nw, w)
          end if
          a(n - 1) = a(0) - a(1)
          a(1) = a(0) + a(1)
          do j = m - 2, 2, -2
              a(2 * j + 1) = a(j) + a(j + 1)
              a(2 * j - 1) = a(j) - a(j + 1)
          end do
          l = 2
          m = mh
          do while (m .ge. 2)
              call dctsub(m, t, nc, w(nw))
              if (m .gt. 4) then
                  call cftfsub(m, t, ip, nw, w)
                  call rftfsub(m, t, nc, w(nw))
              else if (m .eq. 4) then
                  call cftfsub(m, t, ip, nw, w)
              end if
              a(n - l) = t(0) - t(1)
              a(l) = t(0) + t(1)
              k = 0
              do j = 2, m - 2, 2
                  k = k + 4 * l
                  a(k - l) = t(j) - t(j + 1)
                  a(k + l) = t(j) + t(j + 1)
              end do
              l = 2 * l
              mh = m / 2
              do j = 0, mh - 1
                  k = m - j
                  t(j) = t(m + k) - t(m + j)
                  t(k) = t(m + k) + t(m + j)
              end do
              t(mh) = t(m + mh)
              m = mh
          end do
          a(l) = t(0)
          a(n) = t(2) - t(1)
          a(0) = t(2) + t(1)
      else
          a(1) = a(0)
          a(2) = t(0)
          a(0) = t(1)
      end if
      end
!
      subroutine dfst(n, a, t, ip, w)
      integer n, ip(0 : *), j, k, l, m, mh, nw, nc
      real*8 a(0 : n - 1), t(0 : n / 2 - 1), w(0 : *), xr, xi, yr, yi
      nw = ip(0)
      if (n .gt. 8 * nw) then
          nw = n / 8
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 2 * nc) then
          nc = n / 2
          call makect(nc, ip, w(nw))
      end if
      if (n .gt. 2) then
          m = n / 2
          mh = m / 2
          do j = 1, mh - 1
              k = m - j
              xr = a(j) + a(n - j)
              xi = a(j) - a(n - j)
              yr = a(k) + a(n - k)
              yi = a(k) - a(n - k)
              a(j) = xr
              a(k) = yr
              t(j) = xi + yi
              t(k) = xi - yi
          end do
          t(0) = a(mh) - a(n - mh)
          a(mh) = a(mh) + a(n - mh)
          a(0) = a(m)
          call dstsub(m, a, nc, w(nw))
          if (m .gt. 4) then
              call cftfsub(m, a, ip, nw, w)
              call rftfsub(m, a, nc, w(nw))
          else if (m .eq. 4) then
              call cftfsub(m, a, ip, nw, w)
          end if
          a(n - 1) = a(1) - a(0)
          a(1) = a(0) + a(1)
          do j = m - 2, 2, -2
              a(2 * j + 1) = a(j) - a(j + 1)
              a(2 * j - 1) = -a(j) - a(j + 1)
          end do
          l = 2
          m = mh
          do while (m .ge. 2)
              call dstsub(m, t, nc, w(nw))
              if (m .gt. 4) then
                  call cftfsub(m, t, ip, nw, w)
                  call rftfsub(m, t, nc, w(nw))
              else if (m .eq. 4) then
                  call cftfsub(m, t, ip, nw, w)
              end if
              a(n - l) = t(1) - t(0)
              a(l) = t(0) + t(1)
              k = 0
              do j = 2, m - 2, 2
                  k = k + 4 * l
                  a(k - l) = -t(j) - t(j + 1)
                  a(k + l) = t(j) - t(j + 1)
              end do
              l = 2 * l
              mh = m / 2
              do j = 1, mh - 1
                  k = m - j
                  t(j) = t(m + k) + t(m + j)
                  t(k) = t(m + k) - t(m + j)
              end do
              t(0) = t(m + mh)
              m = mh
          end do
          a(l) = t(0)
      end if
      a(0) = 0
      end
!
! -------- initializing routines --------
!
      subroutine makewt(nw, ip, w)
      integer nw, ip(0 : *), j, nwh, nw0, nw1
      real*8 w(0 : nw - 1), delta, wn4r, wk1r, wk1i, wk3r, wk3i
      ip(0) = nw
      ip(1) = 1
      if (nw .gt. 2) then
          nwh = nw / 2
          delta = atan(1.0d0) / nwh
          wn4r = cos(delta * nwh)
          w(0) = 1
          w(1) = wn4r
          if (nwh .eq. 4) then
              w(2) = cos(delta * 2)
              w(3) = sin(delta * 2)
          else if (nwh .gt. 4) then
              call makeipt(nw, ip)
              w(2) = 0.5d0 / cos(delta * 2)
              w(3) = 0.5d0 / cos(delta * 6)
              do j = 4, nwh - 4, 4
                  w(j) = cos(delta * j)
                  w(j + 1) = sin(delta * j)
                  w(j + 2) = cos(3 * delta * j)
                  w(j + 3) = -sin(3 * delta * j)
              end do
          end if
          nw0 = 0
          do while (nwh .gt. 2)
              nw1 = nw0 + nwh
              nwh = nwh / 2
              w(nw1) = 1
              w(nw1 + 1) = wn4r
              if (nwh .eq. 4) then
                  wk1r = w(nw0 + 4)
                  wk1i = w(nw0 + 5)
                  w(nw1 + 2) = wk1r
                  w(nw1 + 3) = wk1i
              else if (nwh .gt. 4) then
                  wk1r = w(nw0 + 4)
                  wk3r = w(nw0 + 6)
                  w(nw1 + 2) = 0.5d0 / wk1r
                  w(nw1 + 3) = 0.5d0 / wk3r
                  do j = 4, nwh - 4, 4
                      wk1r = w(nw0 + 2 * j)
                      wk1i = w(nw0 + 2 * j + 1)
                      wk3r = w(nw0 + 2 * j + 2)
                      wk3i = w(nw0 + 2 * j + 3)
                      w(nw1 + j) = wk1r
                      w(nw1 + j + 1) = wk1i
                      w(nw1 + j + 2) = wk3r
                      w(nw1 + j + 3) = wk3i
                  end do
              end if
              nw0 = nw1
          end do
      end if
      end
!
      subroutine makeipt(nw, ip)
      integer nw, ip(0 : *), j, l, m, m2, p, q
      ip(2) = 0
      ip(3) = 16
      m = 2
      l = nw
      do while (l .gt. 32)
          m2 = 2 * m
          q = 8 * m2
          do j = m, m2 - 1
              p = 4 * ip(j)
              ip(m + j) = p
              ip(m2 + j) = p + q
          end do
          m = m2
          l = l / 4
      end do
      end
!
      subroutine makect(nc, ip, c)
      integer nc, ip(0 : *), j, nch
      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) = cos(delta * nch)
          c(nch) = 0.5d0 * c(0)
          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 cftfsub(n, a, ip, nw, w)
      integer n, ip(0 : *), nw
      real*8 a(0 : n - 1), w(0 : nw - 1)
      if (n .gt. 8) then
          if (n .gt. 32) then
              call cftf1st(n, a, w(nw - n / 4))
              if (n .gt. 512) then
                  call cftrec4(n, a, nw, w)
              else if (n .gt. 128) then
                  call cftleaf(n, 1, a, nw, w)
              else
                  call cftfx41(n, a, nw, w)
              end if
              call bitrv2(n, ip, a)
          else if (n .eq. 32) then
              call cftf161(a, w(nw - 8))
              call bitrv216(a)
          else
              call cftf081(a, w)
              call bitrv208(a)
          end if
      else if (n .eq. 8) then
          call cftf040(a)
      else if (n .eq. 4) then
          call cftx020(a)
      end if
      end
!
      subroutine cftbsub(n, a, ip, nw, w)
      integer n, ip(0 : *), nw
      real*8 a(0 : n - 1), w(0 : nw - 1)
      if (n .gt. 8) then
          if (n .gt. 32) then
              call cftb1st(n, a, w(nw - n / 4))
              if (n .gt. 512) then
                  call cftrec4(n, a, nw, w)
              else if (n .gt. 128) then
                  call cftleaf(n, 1, a, nw, w)
              else
                  call cftfx41(n, a, nw, w)
              end if
              call bitrv2conj(n, ip, a)
          else if (n .eq. 32) then
              call cftf161(a, w(nw - 8))
              call bitrv216neg(a)
          else
              call cftf081(a, w)
              call bitrv208neg(a)
          end if
      else if (n .eq. 8) then
          call cftb040(a)
      else if (n .eq. 4) then
          call cftx020(a)
      end if
      end
!
      subroutine bitrv2(n, ip, a)
      integer n, ip(0 : *), j, j1, k, k1, l, m, nh, nm
      real*8 a(0 : n - 1), xr, xi, yr, yi
      m = 1
      l = n / 4
      do while (l .gt. 8)
          m = m * 2
          l = l / 4
      end do
      nh = n / 2
      nm = 4 * m
      if (l .eq. 8) then
          do k = 0, m - 1
              do j = 0, k - 1
                  j1 = 4 * j + 2 * ip(m + k)
                  k1 = 4 * k + 2 * ip(m + j)
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nm
                  k1 = k1 + 2 * nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nm
                  k1 = k1 - nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nm
                  k1 = k1 + 2 * nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nh
                  k1 = k1 + 2
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 - nm
                  k1 = k1 - 2 * nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 - nm
                  k1 = k1 + nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 - nm
                  k1 = k1 - 2 * nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + 2
                  k1 = k1 + nh
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nm
                  k1 = k1 + 2 * nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nm
                  k1 = k1 - nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 + nm
                  k1 = k1 + 2 * nm
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 - nh
                  k1 = k1 - 2
                  xr = a(j1)
                  xi = a(j1 + 1)
                  yr = a(k1)
                  yi = a(k1 + 1)
                  a(j1) = yr
                  a(j1 + 1) = yi
                  a(k1) = xr
                  a(k1 + 1) = xi
                  j1 = j1 - nm
                  k1 = k1 - 2 * nm
                  xr = a(j1)