> From arpa!SIERRA!BUNEMAN Wed Jul  6 21:36:09 PDT 1988
> Received: by sierra.STANFORD.EDU (3.2/4.7); Wed, 6 Jul 88 21:36:10 PDT
> Date: Wed 6 Jul 88 21:36:09-PDT
> From: Oscar Buneman <BUNEMAN@SIERRA>
> Subject: Fast Hartley Transform
> To: ehg@research.att.com
> Message-Id: <584253369.0.BUNEMAN@SIERRA>
> Mail-System-Version: <SUN-MM(219)+TOPSLIB(128)@SIERRA>
> 
>  6-Jul-88 20:59:28-PDT,4375;000000000000
> Received: from lindy.Stanford.EDU by sierra.STANFORD.EDU (3.2/4.7); Wed, 6 Jul 88 20:59:25 PDT
> From: D5.N60@Forsythe.Stanford.EDU
> Received: by lindy.Stanford.EDU (4.0/4.7); Wed, 6 Jul 88 20:59:04 PDT
> Date: Wed,  6 Jul 88 21:01:02 PDT
> To: buneman@sierra
> 
> 
$NOFLOATCALLS
$STORAGE:2
      SUBROUTINE FHT(F,N)
C
C   The developers of this  program  and  the  underlying  algorithms  are  R.
C Bracewell and O. Buneman of Stanford University. Users  of  the  program  are
C advised to take note of the following:
C
C   The Board of Trustees of the Leland Stanford Junior University  (Stanford)
C is the owner of the patent (U.S. Patent #  4,646,256)  rights  and  copyright
C rights in this Software. Stanford hereby grants to you, the possessor of this
C copy, the right to copy, adapt and use the Software and derivative works  for
C non- Software or  derivative  work  with  all  copyright,  patent  and  other
C proprietary rights notices included on this Software and so long as you grant
C back to Stanford the right to copy, adapt, and use such derivative work.  All
C derivative works must also contain a notice that it is a  derivative  of  the
C Software obtained from Stanford. You  are  not  permitted  to  use,  copy  or
C distribute the Software or derivative works for payment or for any commercial
C use without Stanford's  prior  written  consent.  A  license  for  commercial
C purposes (use or sale) is available from Stanford. Please contact:
C
C         Office of Technology Licensing
C         Stanford University
C         350 Cambridge Avenue, Suite 250
C         Palo Alto, CA 94306
C         (415) 723-0651
C
C   Copyright c 1987 The Board of  Trustees  of  the  Leland  Stanford  Junior
C University. U.S. Patent # 4,646,256. All rights reserved.
C
      LOGICAL LAMBDA
      DIMENSION F(N),HALSEC(10),COSTAB(12),SINTAB(12),LAMBDA(12)
      DATA HALSEC/.54119610,.50979558,.50241929,.50060300,
     &.50015064,.50003765,.50000941,.50000235,.50000059,.50000015/,
     &COSTAB/0.,.70710678,.92387953,.98078528,.99518473,.99879546,
     &.99969882,.99992470,.99998118,.99999529,.99999882,.99999971/,
     &SINTAB/1.,.70710678,.38268343,.19509032,.09801714,.04906767,
     &.02454123,.01227154,.00613588,.00306796,.00153398,.00076699/
     &LAMBDA /.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,
     &.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE./
C ..LEVEL ZERO:...............................
      LEVEL=0
      Q=1./SQRT(FLOAT(N))
      L=1
      M=N/2
      MORL=M
      IU=1
      U=F(IU)
      GO TO 2
1     U=F(IU)
      F(IX)=X
2     IX=IU+M
      X=F(IX)
      IX=IX-MORL
      MORL=M+L-MORL
      F(IX)=Q*(U+X)
      X=Q*(U-X)
      IX=IX+L
      IU=IU+L
      IF(IU.LE.M)GO TO 1
      F(IX)=X
C ..HIGHER LEVELS:............................
3     LEVEL=LEVEL+1
      IU=1
      IV=L+1
      L=L+L
      IF(L.EQ.N) RETURN
      M=MAX0(L,M/2)
      S=0.
      C=1.
      Q=-1.
      MORL=M
C .."LAMBDA" AND "NU" LOOPS:..................
4     LIMIT=IU+M
      U=F(IU)
      V=F(IV)
      GO TO 6
C .."MU" LOOP:................................
5         U=F(IU)
          V=F(IV)
          F(IX)=X
          F(IY)=Y
6         IX=IU+M
          IY=IV+M
          X=C*F(IX)+S*F(IY)
          Y=S*F(IX)-Q*F(IY)
          IX=IX-MORL
          IY=IY-MORL
          MORL=M+L-MORL
          F(IX)=U+X
          F(IY)=V+Y
          X=U-X
          Y=V-Y
          IX=IX+L
          IY=IY+L
          IU=IU+L
          IV=IV+L
          IF(IU.LT.LIMIT) GO TO 5
C ..END OF "MU" LOOP..........................
        F(IX)=X
        F(IY)=Y
        IU=IU+M
        IV=IV+M
        IF(IU.LT.N) GO TO 4
C END OF "NU" LOOP, RESTORE & AUGMENT ANGLE:..
      IU=IU-N+1
      IV=L-IU+2
C ..INCREMENT LAMBDA BY BINARY LOGIC:.........
      K=LEVEL+1
7       K=K-1
        LAMBDA(K)=.NOT.LAMBDA(K)
        IF(.NOT.LAMBDA(K)) GO TO 7
C ..CHECK IF PI/2 HAS BEEN REACHED:...........
      IF(K.EQ.1) GO TO 9
      S=SINTAB(K)
      C=COSTAB(K)
      Q=C
C ..REPLACE TABLE ENTRIES BY NEXT NEEDED:.....
      IF(K.EQ.2) GO TO 4
      K2=K-1
8       K2=K2-1
        IF(LAMBDA(K2)) GO TO 8
      SINTAB(K)=HALSEC(K-2)*(SINTAB(K-1)+SINTAB(K2))
      COSTAB(K)=HALSEC(K-2)*(COSTAB(K-1)+COSTAB(K2))
      GO TO 4
C ..END OF "LAMBDA" LOOP......................
9     LAMBDA(1)=.FALSE.
      IF(LEVEL.LT.3) GO TO 3
C ..RESTORE TRIG TABLES: SWAP SINES & COSINES:
        DO 10 K=3,LEVEL
        Q=COSTAB(K)
        COSTAB(K)=SINTAB(K)
10      SINTAB(K)=Q
      GO TO 3
      END

.