TWO-DIMENSIONAL FT

The program fth() is set up so that the vectors transformed can be either rows or columns of a two-dimensional array. To see how this works, recall that in Fortran a matrix allocated as (n1,n2) can be subscripted as a matrix (i1,i2) or as a long vector (i1 + n1*(i2-1),1), and call sub(x(i1,i2)) passes the subroutine a pointer to the (i1,i2) element. To transform an entire axis, the subroutines ft1axis() and ft2axis() are given. For a two-dimensional FT, we simply call both ft1axis() and ft2axis() in either order. 

# 1D Fourier transform on a 2D data set along the 1-axis
#
subroutine  ft1axis( adj, sign1, n1,n2, cx)
integer i2,          adj,        n1,n2
complex cx(n1,n2)
real sign1
do i2= 1, n2
                call fth( adj, sign1, 1,n1, cx(1,i2))
return; end

 

# 1D Fourier transform on a 2D data set along the 2-axis
#
subroutine  ft2axis( adj, sign2, n1,n2, cx)
integer i1,          adj,        n1,n2
complex cx(n1,n2)
real    sign2
do i1= 1, n1
                call fth( adj, sign2, n1,n2, cx(i1,1))
return; end

I confess that there are faster ways to do things than those I have shown you above. When we are doing many FTs, for example, the overhead calculations done the first time should be saved for use on subsequent FTs, as in the subroutine rocca() included in IEI. Further, manufacturers of computers for heavy numerical use generally design special FT codes for their architecture. Although the basic fast FT used here ingeniously stores its output on top of its input, that feature is not compatible with vectorizing architectures.