#<
# Spectrum
#
# DESCRIPTION
# Obtain averaged amplitude spectra
#
# USAGE
#  Spectrum j2=all <  in.H > out.H
#
# INPUT PARAMETERS
#  j2  - integer [all] number of traces to average at a time
#  n1,n2,n3 - integer  seplib params  (n1<8192)
#  d1,d2,d3 - real seplib params
#  label1 - char*  first axis label 
#  esize -  integer [4] real input data
#                   [8] complex input data	
#
# OUTPUT PARAMETER
# n1,n2    -  integer cube dimension
#		         (data length are padded to a power of 2)
# label1   - char*  "cycles per "//label1
#
# CATEGORY
# Seis:Filter
#
# COMPILE LEVEL
# DISTR
#>
#%
# Author - Peter Mora
# EDIT HISTORY
#	8-25-85	jon	putch o1=0. title="amplitude spectrum"
#      9-22-87       clem   maxnt to 8192, as in cfft1
#      1-26-88       john   handle complex input data
#	1-8-92		jon	dump nt,nx anachronisms
#
# ----------
#  Keyword:  average amplitude spectrum
# ----------
   integer input,output,fetch,infd,outfd,ier,iexpon2
   integer sreed, srite
   integer nt,nx,np,it,n2,n1,i2out,n2out,j2,n,i3,n3,esize
   real data(8192),spec(8192/2+1),dt,df
   complex cdata(8192)
   character label1*80,label*80

   call initpar()
   call doc(SOURCE)
#
# check to see if the data is real or complex and act appropriately
#   (assume real)
   esize=4
   ier=fetch('esize','i',esize)
   if(fetch('n1','i',nt)==0)	call erexit('n1 is required')
   if(nt>8192)			call erexit('n1 too big')
   n=iexpon2(nt)
   if(nt!=n) {
      call putlin('WARNING: data length will be padded to a power of 2')
   }

   nx=1				;ier=fetch('n2','i',nx)
   np=1				;ier=fetch('n3','i',np)
   j2=nx*np			;ier=fetch('j2','i',j2)
   dt=1				;ier=fetch('d1','r',dt)
   df=1/(n*dt)

   label1=' '			;ier=fetch('label1','s',label)
   label1='cycles per '//label

   n2out=(nx*np-1)/j2+1 ; n1=n/2+1
   call putch('n1','i',n1) ;      call putch('n2','i',n2out) 
   n3=1;                          call putch('n3','i',n3)
   call putch('j2 ','i',j2)
   call putch('d1','r',df)
   call putch('o1','r',0.0)
   call putch('label1','s',label1)
   call putch('title','s','Amplitude Spectrum')
   call hclose()

#  loop over all traces

   do i2out=1,n2out {
      do i1=1,n1 { spec(i1)=0. }
      do i2=1,j2 {
        if(esize==4){
         ier=sreed('in',data,4*nt)
	  do it=1,nt { cdata(it)=data(it) } 
                    }
        if(esize==8){
         ier=sreed('in',cdata,8*nt)
                    }
          do it=nt+1,n { cdata(it)=0. }
	 call cfft1(cdata,n,1)
	 do i1=1,n1 { spec(i1)=spec(i1)+cabs(cdata(i1)) }
      }
      do i1=1,n1 { spec(i1)=spec(i1)/j2 }
      ier=srite('out',spec,4*n1)
   }
   end

   function iexpon2(n)
   integer*4 iexpon2
#  compute required length of fourier series for n points such that
#  log2(iexpon2) = nearest integer greater than or equal to n.

   integer n,m,i
   real alog2
   iexpon2=n
   m=alog2(float(n))
   i=2**m
   if(i!=n) iexpon2=2*i
   return
   end


   function alog2(x)

#  log base 2

   real alog2,x
   alog2=alog(x)/alog(2.)
   return
   end