Fourier.frink

View or download Fourier.frink in plain text format


// Routines for Fourier transforms in Frink.
//
// Since different fields of mathematics and engineering use different
// conventions for the Fourier transform, all of these functions allow you
// to (optionally) specify the scaling factor and sign convention.
//
// The Fourier transform implemented by all of these functions can be
// described, in most general terms, as giving element j of the FFT
// of the sequence:
//
// x0, x1, x2, ... , x_n-1
//
// (Where j = 0, 1, ... n-1) as:
//
//                              n - 1
//                               ---                    __
//                    1          \          direction 2 || i j k / n
//    FFT(j) = ----------------   >   x  * e
//              (1-divFactor)/2  /     k
//             n                 ---
//                              k = 0
//
// (Whew!)  Thanks to the excellent program "JavE" for drawing this
// equation.  http://www.jave.de/ 
//
// 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)
//

// Performs the Discrete Fourier Transform of an array of values.
// This is a na�ve implementation and is O(n^2).
// See the excellent exposition at:
//   http://www.spd.eee.strath.ac.uk/~interact/fourier/dft/dftmaths.html
//
// Note that this algorithm is now built into Frink!
DFT[values, divFactor=-1, direction=1] :=
{
   len = length[values]

//   println["Length is $len"]

   results = new array[len]

   mulFactor = 1/len^((1-divFactor)/2)

   s = (direction 2 pi i)/len

   for k = 0 to len-1
   {
      ss = s*k
      sum = 0 * values@0        // This preserves units
      for n = 0 to len-1
         sum = sum + values@n * e^(ss*n)

      results@k = sum*mulFactor  // Optional scaling factor
   }

   return results
}


// Produces the inverse of the DFT given by the DFT function.  In fact, it just
// calls the DFT function with appropriately-reversed parameters.
//
// If you specified the optional second or third arguments for the DFT
// 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.
//
// See the top of this file for information on optional parameters.
InverseDFT[x, divFactor=-1, direction = 1] := DFT[x, -divFactor, -direction]


// Fast Fourier Transform
//
// (Cooley-Tukey decimation-in-time)
// Crandall and Pomerance algorithm 9.5.5
//
// The first argument is an array containing a list of (real or complex) data.
// This algorithm pads the input data out to the next largest power of 2
// (using zeros,) but does not modify the original.
//
// See the top of this file for information on optional parameters.
FFT[x, divFactor=-1, direction=1] :=
{
   x = padAndScramble[x]
   n = length[x]

   g = e^(direction 2 pi i/n)

   m = 1
   while m<n
   {
      mm = n/(2 m)
      for j = 0 to m-1
      {
         a = g^(j mm)
         ii = j
         while ii<n
         {
            i2 = ii+m
            xi = x@ii
            xim = x@i2
            axim = a xim
            x@ii = xi + axim
            x@i2 = xi - axim
            ii = ii + 2 m
         }
      }
      m = 2 m
   }

   // Optional scaling
   mulFactor = 1/n^((1-divFactor)/2)

   if (mulFactor != 1)
   {
      for ii=0 to n-1
         x@ii = x@ii * mulFactor
   }

   return x
}

// 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.  See the top of the file for more information about these
// parameters.
//
// See the top of this file for information on optional parameters.
InverseFFT[x, divFactor=-1, direction = 1] := FFT[x, -divFactor, -direction]


// Returns a copy of the array x scrambled for use in a FFT and padded out
// to the next largest power of 2.
//
// You don't need to call this; it's done by the FFT function.
//
// Modified from Crandall and Pomerance, algorithm 9.5.5 part 3
padAndScramble[x] :=
{
   n = length[x]
   newSize = 2^(bitLength[n]-1)
   if newSize < n
      newSize = newSize * 2

//   println["New size is $newSize"]
   newArray = new array[newSize]

   zeroUnit = 0 x@0             // Maintain units

   for ii = 0 to n-1
      newArray@ii = x@ii
   
   for ii = n to newSize-1      // Do padding
      newArray@ii = zeroUnit

   j = 0
   for ii = 0 to newSize-2
   {
      if ii<j
      {
         temp = newArray@ii
         newArray@ii = newArray@j
         newArray@j = temp
      }
      k = newSize div 2
      while k<=j
      {
         j = j - k
         k = k div 2
      }

      j = j+k
   }

   return newArray
}


View or download Fourier.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17646 days, 9 hours, 38 minutes ago.