c234567890
c This is Fourier analysis and synthesis. 
c It is entirely real; to show how it works I shall start
c with a time function DATA. Every odd sample is the input time
c function and every even sample is 0.0 because the time function is
c real. Since this fast fourier works on the principle of symmetry
c (see MJK's work "computed offshore seismograms") you must consider
c always 2**N points where N is .ge.2. and.le.15 (say).
c use SIGN=-1. to transform from time to frequency
c use SIGN= 1. to inverse transform from say freq. to time
c EXAMPLE STARTS HERE
C234567
      DIMENSION DATA(2048)
      dimension ampli(1024), phase(1024), copha(1024)
      DO 1 I=1,2048
1     DATA(I)=0.0
      DATA(1)=1.0
      DATA(3)= -2.75
      DATA(5)=9.11
      DATA(7)=-1.11
      DATA(9)= 5.33
      DATA(11)=-8.00
      DATA(13)= 4.00
      DATA(15)=-5.0
      data(17)=-3.33
      data(19)=-5.03
      data(21)= 5.03
      data(23)=10.00
      data(25)=10.00
      data(27)=19.00
      data(29)=-7.77
      data(31)= 0.0
      data(33)=1.0
      data(35)=2.0
      data(37)=-7.99
      data(39)=2.0
      data(41)=-3.0
      data(43)=2.0
      data(45)= -3.999
      data(47)=2.0
      data(49)=3.0
      data(51)=3.0
      data(53)=-3.0
      data(55)=3.7
      data(57)=3.9
      data(59)=-.90
      data(61)=90.0
      data(63)=4.0
      N=5
      ND=2**N
      LDATA=2*ND
      WRITE(*,2)(DATA(2*I-1),I=1,ND)
2     FORMAT(1X,F15.3)
3     FORMAT(9X,F15.3)
      WRITE(*,3)(DATA(2*I),I=1,ND)
      WRITE(*,4)(DATA(I),I=1,LDATA)
4     FORMAT(20X,F15.3)
      sign=-1.0
      call wfft(sign,n,data)
      write(*,5)(data(2*i-1),i=1,nd)
      write(*,6)(data(2*i  ),i=1,nd)
 5    format(2x, f10.3)
 6    format(6x, f10.3)
      sign=1.0
      call wfft(sign,n,data)
      write(*,7)(data(2*i-1),i=1,nd)
      write(*,8)(data(2*i  ),i=1,nd)
 7    format(11x, f10.3)
 8    format(21x, f10.3)
      sign=-1.0
      call wfft(sign,n,data)
      do 90 i=1,nd
      ampli(i)=sqrt(data(2*i-1)**2 + data(2*i)**2)
      phase(i)=atan2(data(2*i),data(2*i-1))
 90   continue
c REMEMBER ATAN2 AND NOT ATAN. 
c ATAN2 is same as subroutine BASTAN.   
      write(*,93)(ampli(i),i=1,nd)
      write(*,70)(phase(i),i=1,nd)
 70   format(40x,f10.3)
 93   format(50x,f10.3)
      do 55 i=1,nd
      data(2*i-1)= ampli(i)*cos(phase(i))
      data(2*i )= ampli(i)*sin(phase(i))
 55   continue
      write(*,17)(data(2*i-1),i=1,nd)
      write(*,18)(data(2*i  ),i=1,nd)
 17   format(30x,f10.3)
 18   format(50x,f10.3)
      sign=1.0
      call wfft(sign,n,data)
      write(*,27)(data(2*i-1),i=1,nd)
      write(*,28)(data(2*i  ),i=1,nd)
27    format(9x,f10.3)
 28   format(3x,f10.3)
      call conpha(nd,phase)
      do 56 i=1,nd
      data(2*i-1)=ampli(i)*cos(phase(i))
      data(2*i  )=ampli(i)*sin(phase(i))
 56   continue
      write(*,65)(ampli(i),i=1,nd)
      write(*,66)(phase(i),i=1,nd)
 65   format(30x,f10.3)
 66   format(40x,f10.3)
      write(*,45)(data(2*i-1),i=1,nd)
      write(*,46)(data(2*i  ),i=1,nd)
 45   format(1x,f10.3)
 46   format(9x,f10.3)
      sign=1.0
      call wfft(sign,n,data)
      write(*,35)(data(2*i-1),i=1,nd)
      write(*,36)(data(2*i  ),i=1,nd)
 35   format(33x,f10.3)
 36   format(15x,f10.3)
      STOP
      END
      subroutine wfft(sign,n,data)
      dimension data(2),inter(25)
      lx=2**n
      ldata=2*lx
      flx=lx
      q1=lx
      il=lx
      pi2=6.28318530718
      sigh=sign*pi2/q1
      do 10 i=1,n
      il=il/2
 10   inter(i)=il
      nblokk=1
      do 40 layer=1,n
      nblock=nblokk
      nblokk=nblokk+nblokk
      lblock=lx/nblock
      l2blok= lblock+lblock
      nw=0
      do 40 iblock=1,nblock
      fnw=nw
      aa=sigh*fnw
      w1=cos(aa)
      w2=sin(aa)
      lstart=l2blok*(iblock-1)
      do 20 i=2,lblock,2
      j1=i+lstart
      j=j1-1
      jh=j+lblock
      jh1=jh+1
      q1=data(jh)*w1-data(jh1)*w2
      q2=data(jh)*w2+data(jh1)*w1
      data(jh)=data(j)-q1
      data(jh1)=data(j1)-q2
      data(j  )=data(j )+q1
      data(j1 )=data(j1)+q2
20    continue
      do 29 i=2,n
      idum=inter(i)
      il=idum+idum
      if(nw-idum-il*(nw/il))40,30,30
 30   nw=nw-idum
 29   continue
 40   nw=nw+idum
      nw=0
      do 60 j=1,lx
      if(nw-j)48,47,47
 47   jj=nw+nw+1
      j1=jj+1
      jh1=j+j
      jh=jh1-1
      q1=data(jj)
      data(jj)=data(jh)
      data(jh)=q1
      q1=data(j1)
      data(j1)=data(jh1)
      data(jh1)=q1
 48   continue
      do 49 i=1,n
      idum=inter(i)
      il=idum+idum
      if(nw-idum-il*(nw/il))60,50,50
 50   nw=nw-idum
 49   continue
 60   nw=nw+idum
      if(sign)7,7,11
 7    do 5 k=1,ldata
 5    data(k)=data(k)/flx
 11   mem=0
      return
      end
      subroutine conpha(last,data)
C CONPHA TAKES THE PHASE WHICH VARIED BETWEEN -3.14159 AND +3.14159
C AND WINDS IT UP UPWARDS SO THAT THE RESULTING OUTPUT PHASE VARIES
C BETWEEN -3.14159 AND SAY + BILLION.
      dimension data(last)
      pi2=6.28318530718
      pi=pi2/2.
      pj=0.0
      do 40 i=2, last
      temp=data(i)+pj-data(i-1)
      if(abs(temp)-pi)40,40,10
 10   if(temp)20,40,30
 20   pj=pj +pi2
      go to 40
 30   pj=pj-pi2
 40   data(i)=data(i)+pj
      return
      end
c     subroutine bastan(z1,z2,batan)
c     zer=0.0
c     pi2=6.28318530718
c     py=pi2/2.
c     zha=atan(z2/z1)
c     if(z1.gt.zer)zaz=zha
c     if(z2.lt.zer.and.z1.lt.zer)zaz=zha-py
c     if(z2.lt.zer.and.z1.eq.zer)zaz=-py/2.
c     if(z2.eq.zer.and.z1.lt.zer)zaz=-py
c     if(z2.eq.zer.and.z1.eq.zer)zaz=-py/4.0
c     if(z2.gt.zer.and.z1.lt.zer)zaz=zha+py
c     if(z2.gt.zer.and.z1.eq.zer)zaz=py/2.
c     batan=zaz
c     return
c     end