program fft implicit real*8 (a-h,o-z) character*12 if,of character ans,yes dimension data(20000) yes='y' 10 print *,'FFT forward (of real function) <1> or backward <2>?' read(5,*) idir print *,'input file name?' read(5,15) if 15 format(a12) call input(data,dt,ndata,idir,if) if(ndata.lt.2) go to 10 call prep(data,nbit,ndata) call swap(data,nbit) call calc(data,nbit,idir) print *,'output real/imag into separate files?' read(5,55) ans if(ans.eq.yes) then print *,'output file name for real part?' read(5,15) of call output(data,dt,ndata,of,idir,1) print *,'output file name for imaginary part?' read(5,15) of call output(data,dt,ndata,of,idir,2) else print *,'output file name?' read(5,15) of call output(data,dt,ndata,of,idir,3) endif print *,'start over?' read(5,55) ans 55 format(a1) if(ans.eq.yes) go to 10 stop end c c c subroutine input(data,dt,ndata,idir,if) implicit real*8 (a-h,o-z) dimension data(1) character*12 if open(1,file=if) rewind(1) ndata=0 if(idir.eq.1) then do 20 i=1,10000 j=2*(i-1)+1 read(1,*,end=30) t,data(j) data(j+1)=0.d0 t0=t1 t1=t ndata=ndata+1 20 continue else do 50 i=1,10000 j=2*(i-1)+1 read(1,*,end=30) t,data(j) read(1,*,end=30) t,data(j+1) ndata=ndata+1 t0=t1 t1=t 50 continue endif 30 close(1) dt=t1-t0 return end c c c subroutine output(data,dt,ndata,of,idir,mout) implicit real*8 (a-h,o-z) dimension data(1) character*12 of open(2,file=of) rewind(2) ddt=1.d0/(dfloat(ndata)*dt) if(idir.eq.2) then factor=1.d0/dfloat(ndata) n2=ndata else factor=1.d0 n2=ndata/2 endif if(mout.eq.1) then do 40 i=1,n2 t=(i-1)*ddt j=2*(i-1)+1 write(2,45) t,data(j)*factor 40 continue elseif(mout.eq.2) then do 50 i=1,n2 t=(i-1)*ddt j=2*i write(2,45) t,data(j)*factor 50 continue else do 60 i=1,n2 t=(i-1)*ddt j=2*(i-1)+1 write(2,45) t,data(j)*factor write(2,45) t,data(j+1)*factor 60 continue endif 45 format(1x,1p,e12.5,', ',e12.5) close(2) return end c c c subroutine calc(data,nbit,idir) implicit real*8 (a-h,o-z) dimension data(1) n=2**nbit do 10 i=1,nbit ke=2**i ke1=ke/2 ur=1.d0 ui=0.d0 th=3.141592654/dfloat(ke1) if(idir.eq.2) th=-th wr=dcos(th) wi=dsin(th) do 20 j=1,ke1 do 30 k=j,n,ke kp=k+ke1 mp=2*(kp-1)+1 mk=2*(k-1)+1 tr=data(mp)*ur-data(mp+1)*ui ti=data(mp)*ui+data(mp+1)*ur data(mp)=data(mk)-tr data(mp+1)=data(mk+1)-ti data(mk)=data(mk)+tr data(mk+1)=data(mk+1)+ti 30 continue temp=ur*wr-ui*wi ui=ur*wi+ui*wr ur=temp 20 continue 10 continue return end c c c subroutine prep(data,nbit,ndata) implicit real*8 (a-h,o-z) dimension data(1) do 10 i=1,30 n=2**i if(ndata.le.n) then nbit=i do 20 j=ndata+1,n jj=2*(j-1)+1 data(jj)=0.d0 data(jj+1)=0.d0 20 continue go to 30 endif 10 continue 30 ndata=n return end c c c function nrev(i,nbit) n=0 k=i-1 do 10 l=1,nbit n=2*n+mod(k,2) k=int(k/2) 10 continue nrev=n+1 return end c c c subroutine swap(data,nbit) implicit real*8 (a-h,o-z) dimension data(1) n=2**nbit do 10 i=1,n m=nrev(i,nbit) if(m.gt.i) then j=2*(i-1)+1 k=2*(m-1)+1 tempr=data(j) tempi=data(j+1) data(j)=data(k) data(j+1)=data(k+1) data(k)=tempr data(k+1)=tempi endif 10 continue return end c c c