* fft.f
*
      integer nmax,np,nf,nf2,nb,nb2,nb4
      real sec,dt
      real t(72000),at,tt,tn
      real s(16384),frq(32768),dbt
* fourie parameter
      complex ft(32768)
      real ftr(32768),fti(32768)
* power spectra (single)
      real pt(32768),pt1(32768)
* paramters
      dt=0.000001
      nf=2
      nf2=nf/2
*
* data reading
      at=0.0
      nmax=0
      do 100 j=1,36000,1
        read(*,*,end=199) t(j)
        t(j)=t(j)*1000.0
        at=at+t(j)
        nmax=nmax+1
100   continue
199   at=at/real(nmax)
*
* variance
      tt=0.
      do 200 j=1,nmax,1
        t(j)=t(j)-at
        tt=tt+t(j)*t(j)/real(nmax)
200   continue
*
* fourie transform
      np=alog(real(nmax)+0.5)/alog(2.0)
      nb=2**np
      nb2=nb/2
      nb4=nb/4
      do 300 j=1,nb4-1,1
        s(j)=sin(real(j)*6.283185/real(nb))
300   continue
      call fft(t,ftr,fti,nb,nb2,nb4,np,s)
      fti(1)=0.
      do 400 j=1,nb2,1
        ft(j)=ftr(j)+fti(j)*(0.0,1.0)
400   continue
*
* spectrum
      tn=2.0*dt/real(nb-1)
      dbt=1.0/real(nb-1)/dt*real(nf2)
      do 500 j=1,nb2,1
        pt1(j)=ft(j)*conjg(ft(j))*tn
500   continue
      call filter(pt1,pt,nb2,mf,nf,nf2)
      do 510 j=1,nb,1
        frq(j)=dbt*real(j-1)
510   continue
*
* output
      open(30,file="spectr.dat")
      do 600 j=2,mf,1 
*        write(30,1000) frq(j),pt(j)/tt*frq(j)*dt*1000.
        write(30,1000) frq(j),pt(j)
1000    format(f12.4,f12.8)
600   continue
      close(30)
      end
*********************************************
      subroutine filter(ssp,sp,m,mf,nf,nf2)
      real ssp(32768),sp(32768)
      real sum
      integer m,mf,nf,nf2
      sum=0.
      mf=0
      do 100 k=1,m,nf2
	mf=mf+1
	sum=ssp(k)*nf
	if(nf.eq.1) goto 199
	do 150 i=1,nf-1,1
	  k1=k-i
	  k2=k+i
	  if(k1.le.0) k1=-(k1-1)+1
	  if(k2.ge.m) k2=m-(k2-m)
	  sum=sum+(ssp(k1)+ssp(k2))*real(nf-1)
150	continue
199	sp(mf)=sum/real(nf*nf)
100   continue
      end
*
      subroutine fft(x,a,b,n,n2,n4,np,s)
      real x(72000),a(32768),b(32768),s(16384)
      integer n,n2,n4,np
      integer ip,m1,m0,m2
      real s1,s2
      do 100 i=1,n2,1
	a(i)=x(i)+x(i+n2)
	b(i)=x(i)-x(i+n2)
100   continue
      do 200 i=1,n2,1
	x(2*i-1)=a(i)
	x(2*i)=b(i)
200   continue
      ip=n4
      m1=1
      do 300 m=1,np-1
	m2=2*m1
	m0=m1-1
	do 310 i=1,ip
	  i1=m2*(i-1)+1
	  i2=i1+m1
	  i3=i1+n2
	  i4=i3+m1
	  a(i1)=x(i1)+x(i3)
	  b(i1)=x(i1)-x(i3)
	  a(i2)=x(i2)
	  b(i2)=x(i4)
	  if(m.eq.1) goto 310
	  do 320 k=1,m0,1
	    s1=s(n4-ip*k)*x(i3+k)-s(ip*k)*x(i4+k)
	    s2=s(ip*k)*x(i3+k)+s(n4-ip*k)*x(i4+k)
	    i5=i1+k
	    i6=i1-k+m2
	    a(i5)=x(i5)+s1
	    a(i6)=x(i5)-s1
	    b(i5)=x(i2+k)+s2
	    b(i6)=-x(i2+k)+s2
320	  continue
310	continue
	do 350 j=1,ip,1
	  j1=m2*(j-1)
	  j2=2*j1
	  j3=j2+m2
	  do 360 k=1,m2
	    x(j2+k)=a(j1+k)
	    x(j3+k)=b(j1+k)
360	  continue
350	continue
	if(ip.eq.1) goto 999
	m1=m2
	ip=ip/2
300   continue
999   end