fft4ft.f

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

! test of fft4f.f
!
      include 'fft4f.f'
!
      program main
      integer nmax, nmaxsqrt
      parameter (nmax = 8192)
      parameter (nmaxsqrt = 64)
      integer n, ip(0 : nmaxsqrt + 1)
      real*8 a(0 : nmax), w(0 : nmax * 5 / 4 - 1), t(0 : nmax / 2), 
     &    err, errorcheck
!
      write (*, *) 'data length n=? '
      read (*, *) n
      ip(0) = 0
!
!   check of CDFT
      call putdata(0, n - 1, a)
      call cdft(n, 1, a, ip, w)
      call cdft(n, -1, a, ip, w)
      err = errorcheck(0, n - 1, 2.0d0 / n, a)
      write (*, *) 'cdft err= ', err
!
!   check of RDFT
      call putdata(0, n - 1, a)
      call rdft(n, 1, a, ip, w)
      call rdft(n, -1, a, ip, w)
      err = errorcheck(0, n - 1, 2.0d0 / n, a)
      write (*, *) 'rdft err= ', err
!
!   check of DDCT
      call putdata(0, n - 1, a)
      call ddct(n, 1, a, ip, w)
      call ddct(n, -1, a, ip, w)
      a(0) = a(0) * 0.5d0
      err = errorcheck(0, n - 1, 2.0d0 / n, a)
      write (*, *) 'ddct err= ', err
!
!   check of DDST
      call putdata(0, n - 1, a)
      call ddst(n, 1, a, ip, w)
      call ddst(n, -1, a, ip, w)
      a(0) = a(0) * 0.5d0
      err = errorcheck(0, n - 1, 2.0d0 / n, a)
      write (*, *) 'ddst err= ', err
!
!   check of DFCT
      call putdata(0, n, a)
      a(0) = a(0) * 0.5d0
      a(n) = a(n) * 0.5d0
      call dfct(n, a, t, ip, w)
      a(0) = a(0) * 0.5d0
      a(n) = a(n) * 0.5d0
      call dfct(n, a, t, ip, w)
      err = errorcheck(0, n, 2.0d0 / n, a)
      write (*, *) 'dfct err= ', err
!
!   check of DFST
      call putdata(1, n - 1, a)
      call dfst(n, a, t, ip, w)
      call dfst(n, a, t, ip, w)
      err = errorcheck(1, n - 1, 2.0d0 / n, a)
      write (*, *) 'dfst err= ', err
!
      end
!
!
      subroutine putdata(nini, nend, a)
      integer nini, nend, j, seed
      real*8 a(0 : *), drnd
      seed = 0
      do j = nini, nend
          a(j) = drnd(seed)
      end do
      end
!
!
      function errorcheck(nini, nend, scale, a)
      integer nini, nend, j, seed
      real*8 scale, a(0 : *), drnd, err, e, errorcheck
      err = 0
      seed = 0
      do j = nini, nend
          e = drnd(seed) - a(j) * scale
          err = max(err, abs(e))
      end do
      errorcheck = err
      end
!
!
! random number generator, 0 <= drnd < 1
      real*8 function drnd(seed)
      integer seed
      seed = mod(seed * 7141 + 54773, 259200)
      drnd = seed * (1.0d0 / 259200)
      end
!