/* This file contains routines for 2-dimensional Fourier transforms of discrete
   data.

   Note that these algorithms are now built into Frink and this file is no
   longer necessary!  */

// The (optional) second argument divFactor sets the scaling factor for
// the results.  This means that the scaling factor (the whole
// expression to the left of the summation symbol above) becomes:
//
//                          FFT        InverseFFT
//
//    divFactor =  0:     1/sqrt[n]    1/sqrt[n]
//    divFactor = -1:     1/n          1          (default)
//    divFactor =  1:     1            1/n      
//
//
// The (optional) third argument direction sets the sign used in the
// exponent.
//
//                              FFT          |        InverseFFT
//                                           |
//    direction =  1:   e^( 2 pi i j k / n)  |  e^(-2 pi i j k / n)   (default)
//    direction = -1:   e^(-2 pi i j k / n)  |  e^( 2 pi i j k / n)


/* Return the Fourier transform of a 2-D array. */
DFT2D[values, divFactor=-1, direction=1] :=
{
   rows = length[values]
   cols = length[values@0]
   result1 = new array
   
   for rowNum = 0 to rows-1
      result1@rowNum = DFT[values@rowNum, divFactor, direction]

   result1 = transpose[result1]

   result2 = new array
   for colNum = 0 to cols-1
      result2@colNum = DFT[result1@colNum, divFactor, direction]

   return transpose[result2]
}

// Produces the inverse of the FFT given by the FFT function.  In fact, it just
// calls the FFT function with appropriately-reversed parameters.
//
// If you specified the optional second or third arguments for the FFT
// function, you will need to pass in the *same* arguments to this function
// to get the inverse operation.  This function takes care of reversing them
// appropriately.
InverseDFT2D[x, divFactor=-1, direction=1] := DFT2D[x, -divFactor, -direction]


/** Returns the transpose of a 2D matrix. */
transpose[values] :=
{
   rows = length[values]
   cols = length[values@0]
   result = new array[[cols, rows]]
   for rowNum = 0 to rows-1
      for colNum = 0 to cols-1
         result@colNum@rowNum = values@rowNum@colNum

   return result      
}