fft2f.f

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

          wdr = wi * wi
          wdi = wi * wr
          ss = 4 * wdi
          wmr = 1 - 2 * wdr
          wmi = 2 * wdi
          if (wmi .ge. 0) then
              call cdft(n, wmr, wmi, a)
              xi = a(0) - a(1)
              a(0) = a(0) + a(1)
              a(1) = xi
          end if
          do k = n / 2 - 4, 4, -4
              j = n - k
              xr = a(k + 2) - a(j - 2)
              xi = a(k + 3) + a(j - 1)
              yr = wdr * xr - wdi * xi
              yi = wdr * xi + wdi * xr
              a(k + 2) = a(k + 2) - yr
              a(k + 3) = a(k + 3) - yi
              a(j - 2) = a(j - 2) + yr
              a(j - 1) = a(j - 1) - yi
              wkr = wkr + ss * wdi
              wki = wki + ss * (0.5d0 - wdr)
              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
              wdr = wdr + ss * wki
              wdi = wdi + ss * (0.5d0 - wkr)
          end do
          j = n - 2
          xr = a(2) - a(j)
          xi = a(3) + a(j + 1)
          yr = wdr * xr - wdi * xi
          yi = wdr * xi + wdi * xr
          a(2) = a(2) - yr
          a(3) = a(3) - yi
          a(j) = a(j) + yr
          a(j + 1) = a(j + 1) - yi
          if (wmi .lt. 0) then
              a(1) = 0.5d0 * (a(0) - a(1))
              a(0) = a(0) - a(1)
              call cdft(n, wmr, wmi, a)
          end if
      else
          if (wi .lt. 0) then
              a(1) = 0.5d0 * (a(0) - a(1))
              a(0) = a(0) - a(1)
          end if
          if (n .gt. 2) then
              xr = a(0) - a(2)
              xi = a(1) - a(3)
              a(0) = a(0) + a(2)
              a(1) = a(1) + a(3)
              a(2) = xr
              a(3) = xi
          end if
          if (wi .ge. 0) then
              xi = a(0) - a(1)
              a(0) = a(0) + a(1)
              a(1) = xi
          end if
      end if
      end
!
      subroutine ddct(n, wr, wi, a)
      integer n, j, k, m
      real*8 wr, wi, a(0 : n - 1), wkr, wki, wdr, wdi, ss, xr
      if (n .gt. 2) then
          wkr = 0.5d0
          wki = 0.5d0
          wdr = 0.5d0 * (wr - wi)
          wdi = 0.5d0 * (wr + wi)
          ss = 2 * wi
          if (wi .lt. 0) then
              xr = a(n - 1)
              do k = n - 2, 2, -2
                  a(k + 1) = a(k) - a(k - 1)
                  a(k) = a(k) + a(k - 1)
              end do
              a(1) = 2 * xr
              a(0) = 2 * a(0)
              call rdft(n, 1 - ss * wi, ss * wr, a)
              xr = wdr
              wdr = wdi
              wdi = xr
              ss = -ss
          end if
          m = n / 2
          do k = 1, m - 3, 2
              j = n - k
              xr = wdi * a(k) - wdr * a(j)
              a(k) = wdr * a(k) + wdi * a(j)
              a(j) = xr
              wkr = wkr - ss * wdi
              wki = wki + ss * wdr
              xr = wki * a(k + 1) - wkr * a(j - 1)
              a(k + 1) = wkr * a(k + 1) + wki * a(j - 1)
              a(j - 1) = xr
              wdr = wdr - ss * wki
              wdi = wdi + ss * wkr
          end do
          k = m - 1
          j = n - k
          xr = wdi * a(k) - wdr * a(j)
          a(k) = wdr * a(k) + wdi * a(j)
          a(j) = xr
          a(m) = (wki + ss * wdr) * a(m)
          if (wi .ge. 0) then
              call rdft(n, 1 - ss * wi, ss * wr, a)
              xr = a(1)
              do k = 2, n - 2, 2
                  a(k - 1) = a(k) - a(k + 1)
                  a(k) = a(k) + a(k + 1)
              end do
              a(n - 1) = xr
          end if
      else
          if (wi .ge. 0) then
              xr = 0.5d0 * (wr + wi) * a(1)
              a(1) = a(0) - xr
              a(0) = a(0) + xr
          else
              xr = a(0) - a(1)
              a(0) = a(0) + a(1)
              a(1) = 0.5d0 * (wr - wi) * xr
          end if
      end if
      end
!
      subroutine ddst(n, wr, wi, a)
      integer n, j, k, m
      real*8 wr, wi, a(0 : n - 1), wkr, wki, wdr, wdi, ss, xr
      if (n .gt. 2) then
          wkr = 0.5d0
          wki = 0.5d0
          wdr = 0.5d0 * (wr - wi)
          wdi = 0.5d0 * (wr + wi)
          ss = 2 * wi
          if (wi .lt. 0) then
              xr = a(n - 1)
              do k = n - 2, 2, -2
                  a(k + 1) = a(k) + a(k - 1)
                  a(k) = a(k) - a(k - 1)
              end do
              a(1) = -2 * xr
              a(0) = 2 * a(0)
              call rdft(n, 1 - ss * wi, ss * wr, a)
              xr = wdr
              wdr = -wdi
              wdi = xr
              wkr = -wkr
          end if
          m = n / 2
          do k = 1, m - 3, 2
              j = n - k
              xr = wdi * a(j) - wdr * a(k)
              a(k) = wdr * a(j) + wdi * a(k)
              a(j) = xr
              wkr = wkr - ss * wdi
              wki = wki + ss * wdr
              xr = wki * a(j - 1) - wkr * a(k + 1)
              a(k + 1) = wkr * a(j - 1) + wki * a(k + 1)
              a(j - 1) = xr
              wdr = wdr - ss * wki
              wdi = wdi + ss * wkr
          end do
          k = m - 1
          j = n - k
          xr = wdi * a(j) - wdr * a(k)
          a(k) = wdr * a(j) + wdi * a(k)
          a(j) = xr
          a(m) = (wki + ss * wdr) * a(m)
          if (wi .ge. 0) then
              call rdft(n, 1 - ss * wi, ss * wr, a)
              xr = a(1)
              do k = 2, n - 2, 2
                  a(k - 1) = a(k + 1) - a(k)
                  a(k) = a(k) + a(k + 1)
              end do
              a(n - 1) = -xr
          end if
      else
          if (wi .ge. 0) then
              xr = 0.5d0 * (wr + wi) * a(1)
              a(1) = xr - a(0)
              a(0) = a(0) + xr
          else
              xr = a(0) + a(1)
              a(0) = a(0) - a(1)
              a(1) = 0.5d0 * (wr - wi) * xr
          end if
      end if
      end
!
      subroutine bitrv(n, a)
      integer n, j, k, l, m, m2, n1
      real*8 a(0 : n - 1), xr
      if (n .gt. 2) then
          m = n / 4
          m2 = 2 * m
          n1 = n - 1
          k = 0
          do j = 0, m2 - 2, 2
              if (j .lt. k) then
                  xr = a(j)
                  a(j) = a(k)
                  a(k) = xr
              else if (j .gt. k) then
                  xr = a(n1 - j)
                  a(n1 - j) = a(n1 - k)
                  a(n1 - k) = xr
              end if
              xr = a(j + 1)
              a(j + 1) = a(m2 + k)
              a(m2 + k) = xr
              l = m
              do while (k .ge. l)
                  k = k - l
                  l = l / 2
              end do
              k = k + l
          end do
      end if
      end
!
      subroutine dfct(n, wr, wi, a)
      integer n, j, k, m, mh
      real*8 wr, wi, a(0 : n), wmr, wmi, xr, xi, an
      wmr = wr
      wmi = wi
      m = n / 2
      do j = 0, m - 1
          k = n - j
          xr = a(j) + a(k)
          a(j) = a(j) - a(k)
          a(k) = xr
      end do
      an = a(n)
      do while (m .ge. 2)
          call ddct(m, wmr, wmi, a)
          xr = 1 - 2 * wmi * wmi
          wmi = 2 * wmi * wmr
          wmr = xr
          call bitrv(m, a)
          mh = m / 2
          xi = a(m)
          a(m) = a(0)
          a(0) = an - xi
          an = an + xi
          do j = 1, mh - 1
              k = m - j
              xr = a(m + k)
              xi = a(m + j)
              a(m + j) = a(j)
              a(m + k) = a(k)
              a(j) = xr - xi
              a(k) = xr + xi
          end do
          xr = a(mh)
          a(mh) = a(m + mh)
          a(m + mh) = xr
          m = mh
      end do
      xi = a(1)
      a(1) = a(0)
      a(0) = an + xi
      a(n) = an - xi
      call bitrv(n, a)
      end
!
      subroutine dfst(n, wr, wi, a)
      integer n, j, k, m, mh
      real*8 wr, wi, a(0 : n - 1), wmr, wmi, xr, xi
      wmr = wr
      wmi = wi
      m = n / 2
      do j = 1, m - 1
          k = n - j
          xr = a(j) - a(k)
          a(j) = a(j) + a(k)
          a(k) = xr
      end do
      a(0) = a(m)
      do while (m .ge. 2)
          call ddst(m, wmr, wmi, a)
          xr = 1 - 2 * wmi * wmi
          wmi = 2 * wmi * wmr
          wmr = xr
          call bitrv(m, a)
          mh = m / 2
          do j = 1, mh - 1
              k = m - j
              xr = a(m + k)
              xi = a(m + j)
              a(m + j) = a(j)
              a(m + k) = a(k)
              a(j) = xr + xi
              a(k) = xr - xi
          end do
          a(m) = a(0)
          a(0) = a(m + mh)
          a(m + mh) = a(mh)
          m = mh
      end do
      a(1) = a(0)
      a(0) = 0
      call bitrv(n, a)
      end
!