fft4f.f

FFT代原碼為C++需要測過才能用

          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
      m = n / 2
      xr = a(0) + a(n)
      a(0) = a(0) - a(n)
      t(0) = xr - a(m)
      t(m) = xr + a(m)
      if (n .gt. 2) then
          mh = m / 2
          do j = 1, mh - 1
              k = m - j
              xr = a(j) + a(n - j)
              a(j) = a(j) - a(n - j)
              xi = a(k) + a(n - k)
              a(k) = a(k) - a(n - k)
              t(j) = xr - xi
              t(k) = xr + xi
          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) call bitrv2(m, ip(2), a)
          call cftsub(m, a, w)
          if (m .gt. 4) call rftsub(m, a, nc, w(nw))
          xr = a(0) + a(1)
          a(n - 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
          a(1) = xr
          l = 2
          m = mh
          do while (m .ge. 2)
              call dctsub(m, t, nc, w(nw))
              if (m .gt. 4) call bitrv2(m, ip(2), t)
              call cftsub(m, t, w)
              if (m .gt. 4) call rftsub(m, t, nc, w(nw))
              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
      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)
              a(j) = a(j) + a(n - j)
              xi = a(k) - a(n - k)
              a(k) = a(k) + a(n - k)
              t(j) = xr + xi
              t(k) = xr - xi
          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) call bitrv2(m, ip(2), a)
          call cftsub(m, a, w)
          if (m .gt. 4) call rftsub(m, a, nc, w(nw))
          xr = a(0) + a(1)
          a(n - 1) = a(1) - a(0)
          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
          a(1) = xr
          l = 2
          m = mh
          do while (m .ge. 2)
              call dstsub(m, t, nc, w(nw))
              if (m .gt. 4) call bitrv2(m, ip(2), t)
              call cftsub(m, t, w)
              if (m .gt. 4) call rftsub(m, t, nc, w(nw))
              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 : *), 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 cftsub(n, a, w)
      integer n, j, j1, j2, j3, k, k1, ks, l, m
      real*8 a(0 : n - 1), w(0 : *)
      real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      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) + a(j1)
              x0i = a(j + 1) + a(j1 + 1)
              x1r = a(j) - a(j1)
              x1i = a(j + 1) - a(j1 + 1)
              x2r = a(j2) + a(j3)
              x2i = a(j2 + 1) + a(j3 + 1)
              x3r = a(j2) - a(j3)
              x3i = a(j2 + 1) - a(j3 + 1)
              a(j) = x0r + x2r
              a(j + 1) = x0i + x2i
              a(j2) = x0r - x2r
              a(j2 + 1) = x0i - x2i
              a(j1) = x1r - x3i
              a(j1 + 1) = x1i + x3r
              a(j3) = x1r + x3i
              a(j3 + 1) = 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) + a(j1)
                  x0i = a(j + 1) + a(j1 + 1)
                  x1r = a(j) - a(j1)
                  x1i = a(j + 1) - a(j1 + 1)
                  x2r = a(j2) + a(j3)
                  x2i = a(j2 + 1) + a(j3 + 1)
                  x3r = a(j2) - a(j3)
                  x3i = a(j2 + 1) - a(j3 + 1)
                  a(j) = x0r + x2r
                  a(j + 1) = x0i + x2i
                  a(j2) = x2i - x0i
                  a(j2 + 1) = x0r - x2r
                  x0r = x1r - x3i
                  x0i = x1i + x3r
                  a(j1) = wk1r * (x0r - x0i)
                  a(j1 + 1) = wk1r * (x0r + x0i)
                  x0r = x3i + x1r
                  x0i = x3r - x1i
                  a(j3) = wk1r * (x0i - x0r)
                  a(j3 + 1) = 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) + a(j1)
                      x0i = a(j + 1) + a(j1 + 1)
                      x1r = a(j) - a(j1)
                      x1i = a(j + 1) - a(j1 + 1)
                      x2r = a(j2) + a(j3)
                      x2i = a(j2 + 1) + a(j3 + 1)
                      x3r = a(j2) - a(j3)
                      x3i = a(j2 + 1) - a(j3 + 1)
                      a(j) = x0r + x2r
                      a(j + 1) = x0i + x2i
                      x0r = x0r - x2r
                      x0i = x0i - x2i
                      a(j2) = wk2r * x0r - wk2i * x0i
                      a(j2 + 1) = wk2r * x0i + wk2i * x0r
                      x0r = x1r - x3i
                      x0i = x1i + x3r
                      a(j1) = wk1r * x0r - wk1i * x0i
                      a(j1 + 1) = wk1r * x0i + wk1i * x0r
                      x0r = x1r + x3i
                      x0i = x1i - x3r
                      a(j3) = wk3r * x0r - wk3i * x0i
                      a(j3 + 1) = 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) - a(j1)
              x0i = a(j + 1) - a(j1 + 1)
              a(j) = a(j) + a(j1)
              a(j + 1) = a(j + 1) + a(j1 + 1)
              a(j1) = x0r
              a(j1 + 1) = x0i
          end do
      end if
      end
!
      subroutine rftsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi
      ks = 4 * nc / n
      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) - a(j)
          xi = a(k + 1) + a(j + 1)
          yr = wkr * xr - wki * xi
          yi = wkr * xi + wki * xr
          a(k) = a(k) - yr
          a(k + 1) = a(k + 1) - yi
          a(j) = a(j) + yr
          a(j + 1) = a(j + 1) - yi
      end do
      end
!
      subroutine dctsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr
      ks = nc / n
      kk = ks
      m = n / 2
      do k = 1, m - 1
          j = n - k
          wkr = c(kk) - c(nc - kk)
          wki = c(kk) + c(nc - kk)
          kk = kk + ks
          xr = wki * a(k) - wkr * a(j)
          a(k) = wkr * a(k) + wki * a(j)
          a(j) = xr
      end do
      a(m) = 2 * c(kk) * a(m)
      end
!
      subroutine dstsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr
      ks = nc / n
      kk = ks
      m = n / 2
      do k = 1, m - 1
          j = n - k
          wkr = c(kk) - c(nc - kk)
          wki = c(kk) + c(nc - kk)
          kk = kk + ks
          xr = wki * a(j) - wkr * a(k)
          a(j) = wkr * a(j) + wki * a(k)
          a(k) = xr
      end do
      a(m) = 2 * c(kk) * a(m)
      end
!