C
c *** Fourier Transform using Osaki's Textbook(1994) based on single precision
c
      PROGRAM tran1s
C
      PARAMETER(MWV=66000)
c      CHARACTER INFILE*30,head1*10,head2*10,head3*10,head4*10
      CHARACTER INFILE*30
      dimension time(MWV)
      dimension XR1(MWV),XR2(MWV),XR3(MWV)
     *       ,sF1(MWV),sF2(MWV),sF3(MWV),sG(MWV)
      COMPLEX FX1(MWV),FX2(MWV),FX3(MWV)
C
c  Prompt for input RECON'ed file
C
      WRITE(*,1010) 
 1010 FORMAT(1X,/'Enter input file of waveform data')
      READ(*,1020) INFILE
 1020 FORMAT(A30)
C      
c *** Read input data file  ***
C
      OPEN(11,FILE=INFILE)
c        read(11,*)head1,head2,head3,head4
        read(11,*)
        i=1
   10   read(11,*,end=11)time(i),xr1(i),xr2(i),xr3(i)
        i=i+1
        goto 10
   11 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)
      write(*,*)'NT and DT=',NT,DT
C
c *** FOURIER TRANSFORM and Integration or Differenciation ***
C
      NF=2
      DO 200 I=2,16
        NF=2*NF
        IF(NF.GE.NT)GOTO 20
  200 CONTINUE
      WRITE(*,*)' NT>2**16 !!!!!'
      STOP
   20 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF=',IF1,NF
C
      DO 220 I=1,NT
        FX1(I)=CMPLX(XR1(I),0.0)
        FX2(I)=CMPLX(XR2(I),0.0)
        FX3(I)=CMPLX(XR3(I),0.0)
  220 CONTINUE
      DO 240 I=NT+1,NF
        FX1(I)=(0.0,0.0)
        FX2(I)=(0.0,0.0)
        FX3(I)=(0.0,0.0)
  240 CONTINUE
C
      CALL WVFAST(NF,FX1,MWV,-1)
      CALL WVFAST(NF,FX2,MWV,-1)
      CALL WVFAST(NF,FX3,MWV,-1)
      DUR=FLOAT(NF)*DT
      DO 260 I=1,NF
        FX1(I)=FX1(I)*DT
        FX2(I)=FX2(I)*DT
        FX3(I)=FX3(I)*DT
  260 CONTINUE
C
c      FOM=0.0
      DF=1.0/DUR
      DOM=6.283185*DF
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/6.283185
C
c   * Integration or differenciation *
c
      WRITE(*,*)
      WRITE(*,*)' Integration/Differentiation:'
      WRITE(*,*)' NO Change         :                 (Enter=0)'
      WRITE(*,*)' Integration       : Ex.) acc -> vel (Enter=1)'
      WRITE(*,*)' Double Integration: Ex.) acc -> dis (Enter=2)'
      WRITE(*,*)' Differentiation:  : Ex.) Vel -> acc (Enter=3)'
      READ(*,*)NIT
c
      if((NIT.eq.1).or.(NIT.eq.2))then
        write(*,*)
        write(*,*)'Note: Fourier Amp. at 0(Hz) is neglected!!!'
      end if
c
      IF(NIT.EQ.1)THEN
        DO 320 I = 2,NOM
           OMEGA=FLOAT(I-1)*DOM
	   FX1(I) = FX1(I)/cmplx(0.0,omega)
	   FX2(I) = FX2(I)/cmplx(0.0,omega)
	   FX3(I) = FX3(I)/cmplx(0.0,omega)
  320   CONTINUE
      END IF
c
      IF(NIT.EQ.2)THEN
        DO 340 I = 2,NOM
           OMEGA=FLOAT(I-1)*DOM
	   FX1(I) = FX1(I)/(-omega*omega)
	   FX2(I) = FX2(I)/(-omega*omega)
	   FX3(I) = FX3(I)/(-omega*omega)
  340   CONTINUE
      END IF
c
      IF(NIT.EQ.3)THEN
        DO 360 I = 1,NOM
           OMEGA=FLOAT(I-1)*DOM
	   FX1(I) = FX1(I)*cmplx(0.0,omega)
	   FX2(I) = FX2(I)*cmplx(0.0,omega)
	   FX3(I) = FX3(I)*cmplx(0.0,omega)
  360   CONTINUE
      END IF
c
c *** Smoothed Spectra by Parzen's Windows ***
c
      WRITE(*,*)
      WRITE(*,*)'Smoothing Spectra by Parzen Windows'
c
      WRITE(*,*)
      write(*,*)'Delta Freq(Hz)=',DF
      write(*,*)'Enter Band Width(Hz):'
        Band1=1.854305/0.5*DF
        Band2=1.854305*50.0*DF
        write(*,*)'   BAND must be larger than',Band1
        write(*,*)'   BAND must be smaller than',Band2
        write(*,*)'   if Band=0, No smoothing.'
        write(*,*)
      write(*,*)'Enter Band Width(Hz):'
      read(*,*)BAND
      do I=1,NF
          sF1(I)=cabs(FX1(I))
          sF2(I)=cabs(FX2(I))
          sF3(I)=cabs(FX3(I))
      end do
      call SWIN(NOM,sF1,sG,MWV,100,DF,BAND)
      call SWIN(NOM,sF2,sG,MWV,100,DF,BAND)
      call SWIN(NOM,sF3,sG,MWV,100,DF,BAND)
C
c   * Output Original and Smoothed Fourier Spectrum *
c
      write(*,*)
      write(*,*)'Output of Original and Smoothed Fourier Spectrum:'
      WRITE(*,2000)
 2000 FORMAT(/1X,'Enter output file name (csv format)')
      READ(*,2020) INFILE
 2020 FORMAT(A30)
      OPEN(10,FILE=INFILE,STATUS='unknown')
      write(10,*)'freq(Hz),Original X,Smoothed X,Original Y,Smoothed Y,O
     *riginal Z,Smoothed Z'
      DO 400 I = 1,NOM
        fq=FLOAT(I-1)*DF
        write(10,2040)fq,cabs(FX1(I)),sF1(I),cabs(FX2(I)),sF2(I)
     *                  ,cabs(FX3(I)),sF3(I)
  400 CONTINUE
 2040 format(f8.4,6(',',e14.6))
      close(10)
c
        STOP
	END
C
c *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM (Osaki, 1994) ***
C
      SUBROUTINE WVFAST(N,A,MOM,IND)
C
c *** N   : NUMBER OF FOURIER COEFF. OR WAVE ( = 2**I )
c *** A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE ( REAL PART )
c    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
c *** MOM : N IN MAIN DIMMENSION
c         : IF MWV > 2**13(=8192), CAN NOT USE !!!
c *** IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX A(MOM),TEMP,THETA
C
c *** CALUCURATION OF FFT
C
      PAI=3.141593
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(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=CMPLX(0.0,PAI*FLOAT(IND*(K-1))/FLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*EXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
c
c
c *** Soomthed Spectra by Parzens's windows (Osaki, 1994) ***
c
      SUBROUTINE SWIN(NF,F,G,ND,IND,DF,BAND)
c
      dimension F(ND),G(ND),W(101),G1(4497),G2(4497)
c
c   * Initialization *
c
      T=1.0/DF
      if(Band.eq.0.0)goto 120
      UDF=1.854305/Band*DF
      if(UDF.gt.0.5)goto 230
      Lmax=IFIX(2.0/UDF)+1
      if(LMAX.gt.101)goto 240
c
c   * Spectral Window *
c
      W(1)=0.75*UDF
      do 110 L=2,LMAX
        DIF=1.570796*real(L-1)*UDF
        W(L)=W(1)*(sin(DIF)/DIF)**4
  110 continue
c
c   * Conversion from Fourier to Power Spectrum *
c
  120 if(IND.ne.100)goto 140
      G(1)=F(1)**2/T
      do 130 K=2,NF-1
        G(K)=2.0*F(K)**2/T
  130 continue
      G(NF)=F(NF)**2/T
c
c   * Smoothing of Power Spectrum *
c
  140 if(Band.eq.0)goto 210
      LL=Lmax*2-1
      LN=LL-1+NF
      LT=(LL-1)*2+NF
      LE=LT-Lmax+1
      do 150 K=1,LT
        G1(K)=0.0
  150 continue
      do 160 k=1,NF
        G1(LL-1+K)=G(K)
  160 continue
      do 180 K=Lmax,LE
        S=W(1)*G1(K)
        Do 170 L=2,Lmax
          S=S+W(L)*(G1(K-L+1)+G1(K+L-1))
  170   continue
        G2(K)=S
  180 continue
      do 190 L=2,Lmax
        G2(LL+L-1)=G2(LL+L-1)+G2(LL-L+1)
        G2(LN-L+1)=G2(LN-L+1)+G2(LN+L-1)
  190 continue
      do 200 K=1,NF
        G(K)=G2(LL-1+K)
  200 continue
c
c   * Smoothed Fourier Spectrum *
c
  210 F(1)=sqrt(G(1)*T)
      do 220 K=2,NF-1
        F(K)=sqrt(G(K)*T/2.0)
  220 continue
      F(NF)=sqrt(G(NF)*T)
      return
c
  230 write(*,601)
      write(*,*)'Enter 0 for stop'
      read(*,*)Nstop
      if(Nstop.eq.0)stop
  240 write(*,602)
      write(*,*)'Enter 0 for stop'
      read(*,*)Nstop
      if(Nstop.eq.0)stop
c
c   * Format Statement *
c
  601 format('Bandwidth is too narrow!')
  602 format('Bandwidth is too wide!')

      end

