*
* Kiso Enshu 2002, Sample program #6
************************************************************************
      program fft
*
*     Fourier spectrum
*
      implicit none
*
      real        a1
      integer     nmax, nyq, nwrk
      parameter ( a1 = 0.8 )
      parameter ( nmax = 500, nyq = nmax/2, nwrk=2*nmax+15 )
*
      real*4      x( nmax )
      real*4      wrk( nwrk )
      real*4      frq( nyq )
      real*4      pow( nyq )
      real*4      y, w, pi
      integer     i, k, k1, k2, iws
*
      character   cfin*90
      character   cfout1*90
      character   cfout2*90
*
      data        cfin
     &/'/home/hiro/kiso/data/random500.grd'/
      data        cfout1
     &/'/home/hiro/kiso/data/ar500_sp.grd'/
      data        cfout2
     &/'/home/hiro/kiso/data/sp500.grd'/
*
*     open file
      open( 10, file=cfin, form='unformatted',
     &          access='direct', recl=4 )
      open( 20, file=cfout1, form='unformatted',
     &          access='direct', recl=4,
     &          status='unknown'        )
      open( 30, file=cfout2, form='unformatted',
     &          access='direct', recl=4,
     &          status='unknown'        )

      do k = 1, nyq
         frq( k ) = float( k ) / float( nmax ) !! frequency
      enddo

      y = 0.
      do i = 1, nmax
         read( 10, rec=i ) w
         y = a1 * y + w
         x( i ) = y

         write( 20, rec=i ) x( i )
      enddo
      close( 20 )

      call rffti( nmax, wrk )
      call rfftf( nmax, x, wrk )

      do k = 1, nmax
         x( k ) = x( k ) / float( nmax ) !! factorize
      enddo

      i = mod( nmax,2 )
      if( i .eq. 0 ) then       !! even
         do k = 1, nyq-1
            k1 = 2 * k
            k2 = 2 * k + 1
            pow( k ) = x( k1 )**2 + x( k2 )**2
         enddo
         pow( nyq ) = x( 2*nyq )**2
      else                      !! odd
         do k = 1, nyq
            k1 = 2 * k
            k2 = 2 * k + 1
            pow( k ) = x( k1 )**2 + x( k2 )**2
         enddo
      endif

      do k = 1, nyq
         write( 30, rec=k ) pow( k )
      enddo
      close( 30 )
*
*     Dennou
*
      write( *,* ) ' WORKSTATION ID (I)  ?; '
      call sgpwsn
      read( *,* ) iws
      call gropn( iws )
      call grfrm

      call usspnt( nyq, frq, pow )
      call uspfit
      call grstrf

      call ussttl( 'frequency', ' ', 'spectrum density', ' ' )
      call usdaxs

      call uuslnt( 1 )
      call uuslni( 9 )
      call uulin( nyq, frq, pow )

      call grcls

      stop
      end