C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFFT(N,FFT,NDFFT,IND)
C
C *** N     : NUMBER OF FOURIER COEFF. OR WAVE (= 2**I)
C *** FFT(I): FOURIER COEFFICIENTS (N/2+1) OR WAVE (REAL PART)
C    (FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0)
C *** NDFFT : N IN MAIN DIMMENSION
C *** IND   : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 FFT(NDFFT),TEMP,THETA
C
      PAI=4.D0*DATAN(1.D0)
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=FFT(J)
        FFT(J)=FFT(I)
        FFT(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=DCMPLX(0.0D0,PAI*DFLOAT(IND*(K-1))/DFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=FFT(J)*DCEXP(THETA)
          FFT(J)=FFT(I)-TEMP
          FFT(I)=FFT(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
