fft4f2d.f

2维fft程序

                  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