MODULE HartleyFFT_2D
IMPLICIT NONE
PUBLIC  :: fht2d

! From hamill@mailhub.nmg.sms.siemens.com Wed Apr 24 16:21:31 1996
 
! Code converted using TO_F90 by Alan Miller
! Date: 2004-02-04  Time: 13:01:18

! I see from the newsgroups that you are interested in multidimensional
! Hartley transforms.  The attached FORTRAN code works (I think ... but
! it's been a few years since I used it) for 2D and 1D.  For 3D see H.Hao
! and R.N.Bracewell, "A Three-Dimensional DFT Algorithm Using the Fast
! Hartley Transform," Proc IEEE vol. 75 no. 2 (Feb, 87) p.264-266.

! I do not vouch for the correctness of the information or the code, but
! perhaps it is a useful starting point.

! Good luck!

! Jim Hamill

!------------------------------------------------
! hamill@mailhub.nmg.sms.siemens.com
! James Hamill
! Siemens Medical Systems, Nuclear Medicine Group
! 2501 N. Barrington Rd.
! Hoffman Estates, IL  60195
! voice:(847)304-7760       FAX:(847)304-7706
!------------------------------------------------

CONTAINS


!------------------------------------------------------------------------------
!   FHT2D
!        A subroutine for the 2-dimensional Hartley transform of an
!        array  F2(N,N)  where N = 2**P
!        This routine call the 1-d Hartley transform routine FHTFR, for
!        rows and then for columns.
!-------------------------------------------------------------------------------

SUBROUTINE fht2d(nx, f2)
INTEGER, INTENT(IN)   :: nx
REAL, INTENT(IN OUT)  :: f2(0:nx-1,0:nx-1)

REAL     :: a, b, c, d, e, f(0:255)
INTEGER  :: i, j, n, ny, n1, p, u, v

p = INT( LOG(REAL(nx))/LOG(2.) + 0.5)
n  = 2**p
n1 = n - 1
DO j = 0,n1
  CALL fhtfr(p, f2(0,j))
END DO
DO i = 0,n1
  DO j = 0,n1
    f(j)=f2(i,j)
  END DO
  CALL fhtfr(p, f)
  DO j = 0,n1
    f2(i,j)=f(j)
  END DO
END DO

!  UP TO THIS POINT F2(U,V) NOW CONTAIN THE SEPARABLE
!                    CAS(2*pi*u*x)*CAS(2*pi*v*y)   transform

!  The following replacements form the nonseparable 2D-Hartley transform
!    Ref:   RN Bracewell, et al., "Fast Two-dimensional Hartley transform,"
!           Proc IEEE, 74:9,1282 (Sept 1986)

ny = nx
DO v = 1,ny/2-1
  DO u = 1,nx/2-1
    a = f2(u,v)
    b = f2(nx-u,v)
    c = f2(u,ny-v)
    d = f2(nx-u,ny-v)
    e = 0.5*( (a+d)-(b+c) )
    f2(u,v)       = a - e
    f2(nx-u,v)    = b + e
    f2(u,ny-v)    = c + e
    f2(nx-u,ny-v) = d - e
  END DO
END DO

RETURN
END SUBROUTINE fht2d


!-----------------------------------------------------------------------------
!   SUBROUTINE FHTFOR TAKES INPUT F() AND RETURNS THE DHT IN THE SAME F()
!   LENGTH OF F() IS 2**P
!-----------------------------------------------------------------------------
!     THIS IS THE CORRECTED VERSION OF FHTFOR.FOR
!     RN BRACEWELL,  THE HARTLEY TRANSFORM, Oxford Univ. Press (1986)
!     pp. 140-144
!-----------------------------------------------------------------------------

SUBROUTINE fhtfr(p, f)
INTEGER, INTENT(IN)   :: p
REAL, INTENT(IN OUT)  :: f(0:255)

REAL     :: h9, s9(256), t, t9(256)
INTEGER  :: d9, e9, i, j, k9, l9, n, n9, q9, sx9, u9, x9, y9,   &
    a9(0:64), m9(0:10), c9(0:10), v9(0:10)
INTEGER  :: pold = -99
REAL, PARAMETER     :: pi = 3.1415926536

IF (p == 1) THEN
  j    = f(0) + f(1)
  f(1) = f(0) - f(1)
  f(0) = j
  GO TO 9636   !   RETURN
END IF
IF(p == pold) GO TO 9400
n9 = 2**(p-2)
n = 4*n9
c9(5) = n-1
c9(6) = p-1

!   GET POWERS OF 2
i = 1
m9(0) = 1
m9(1) = 2
9204   m9(i+1) = m9(i) + m9(i)
i = i + 1
IF(i <= p) GO TO 9204

!  TRIGONOMETRIC FUNCTION TABLES
IF(n == 4) GO TO 9400
IF(n == 8) THEN
  s9(2) = 1.
  s9(1) = SIN(pi/4.)
  GO TO 9330
END IF
!   GET SINES
s9(n9) = 1.
!   COARSE SEED TABLE FOR SINES
DO i = 1,3
  s9(i*n9/4) = SIN(i*pi/8.)
END DO
!   INITIAL HALF SECANT
h9 = 1./2/COS(pi/16)
!   FILL SINE TABLE
c9(4) = p-4
DO i = 1,p-4
  c9(4) = c9(4)-1
  v9(0) = 0
  DO j = m9(c9(4)), n9-m9(c9(4)), m9(c9(4)+1)
    v9(1) = j + m9(c9(4))
    s9(j) = h9*(s9(v9(1))+v9(0))
    v9(0) = s9(v9(1))
  END DO
!   HALF SECANT RECURSION
  h9 = 1/SQRT(2+1/h9)
END DO

!   GET TANGENTS
9330 c9(0) = n9 - 1
DO i = 1,n9-1
  t9(i) = (1-s9(c9(0)))/s9(i)
  c9(0) = c9(0) - 1
END DO
t9(n9) = 1

!   FAST PERMUTE
!   FOR P = 2,3 PERMUTE DIRECTLY
9400 IF(p == 2) THEN
  v9(9) = f(1)
  f(1)  = f(2)
  f(2)  = v9(9)
  GO TO 9500
END IF
IF(p == 3) THEN
  v9(9) = f(1)
  f(1)  = f(4)
  f(4)  = v9(9)
  v9(9) = f(3)
  f(3)  = f(6)
  f(6)  = v9(9)
  GO TO 9500
END IF
IF(p == pold) GO TO 9420
!   FOR P=4, 5, 6 (Q9=2,3), SKIP STRUCTURE TABLE
q9 = p/2
c9(2) = m9(q9)
q9 = q9 + MOD(p,2)
IF(q9 == 2) THEN
  a9(0) = 0
  a9(1) = 2
  a9(2) = 1
  a9(3) = 3
  GO TO 9420
END IF
IF(q9 == 3) THEN
  a9(0) = 0
  a9(1) = 4
  a9(2) = 2
  a9(3) = 6
  a9(4) = 1
  a9(5) = 5
  a9(6) = 3
  a9(7) = 7
  GO TO 9420
END IF

!   SET UP STRUCTURE TABLE
a9(0) = 0
a9(1) = 1
DO i = 2,q9
  DO j = 0,m9(i-1)-1
    a9(j) = a9(j) + a9(j)
    a9(j+m9(i-1))=a9(j)+1
  END DO
END DO

!   PERMUTE
9420 DO i = 1,c9(2)-1
  v9(4)  = c9(2)*a9(i)
  v9(5)  = i
  v9(6)  = v9(4)
  v9(7)  = f(v9(5))
  f(v9(5)) = f(v9(6))
  f(v9(6)) = v9(7)
  DO j =1,a9(i)-1
    v9(5) = v9(5) + c9(2)
    v9(6) = v9(4) + a9(j)
    v9(7) = f(v9(5))
    f(v9(5)) = f(v9(6))
    f(v9(6)) = v9(7)
  END DO
END DO

!   STAGES 1 & 2
!   GET TWO-ELEMENT DHTs
9500 DO i = 0,n-2,2
  t      = f(i) - f(i+1)
  f(i)   = f(i) + f(i+1)
  f(i+1) = t
END DO

!   GET FOUR-ELEMENT DHTs
DO i = 0,n-4,4
  t      = f(i) - f(i+2)
  f(i)   = f(i) + f(i+2)
  f(i+2) = t
  t      = f(i+1)-f(i+3)
  f(i+1) = f(i+1)+f(i+3)
  f(i+3) = t
END DO

IF(p == 2) GO TO 9636  !  RETURN

!   STAGES 3, 4, ...
u9 = c9(6)
sx9 = 4
DO l9 = 2,c9(6)
  v9(2) = sx9 + sx9
  u9 = u9-1
  v9(3) = m9(u9-1)
  DO q9 = 0,c9(5),v9(2)
    i = q9
    d9 = i+sx9
    t     = f(i) - f(d9)
    f(i)  = f(i) + f(d9)
    f(d9) = t
    k9 = d9-1
    DO j=v9(3),n9,v9(3)
      i  = i  + 1
      d9 = i  + sx9
      e9 = k9 + sx9
      t  = f(d9) + f(e9)*t9(j)
      x9 = f(e9) - t*s9(j)
      y9 = x9*t9(j) + t
      t = f(i) + y9
      f(d9) = f(i)  - y9
      f(e9) = f(k9) + x9
      f(k9) = f(k9) - x9
      f(i)  = t
      k9 = k9 -1
    END DO
  END DO
  sx9= v9(2)
END DO

!   SAVE P SO TABLE GENERATION CAN BE SKIPPED ON NEXT CALL
9636 pold = p
RETURN
END SUBROUTINE fhtfr

END MODULE HartleyFFT_2D