Fourier.frink

Download or view 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.
//
// Note that this algorithm is now built into Frink!
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
}


/** These functions are intended to help you figure out what frequency or
    period each component in the Fourier transform corresponds to. */


/** Calculate the frequency of the specified index given:
    
    [ft, totalPeriod, index]

    where:
      ft: the transformed sequence
      totalPeriod:  the total period of the initial sequence (e.g. 1 year)
      index:  the index for which we want to calculate the index.

    Note that this gives the frequency.  If you want the *angular* frequency
    (omega) such as to pass to a sin or cos function, you need to multiply this
    frequency by 2 pi.
*/

frequencyFromTotalPeriod[ft, totalPeriod, index] :=
{
   if index == 0  // DC component
      return 0 * 1/totalPeriod  // Get the units right?

   return 1/periodFromTotalPeriod[ft, totalPeriod, index]
}


/** Calculate the period of the specified index given:
    
    [ft, totalPeriod, index]

    where:
      ft: the transformed sequence
      totalPeriod:  the total period of the initial sequence (e.g. 1 year)
      index:  the index for which we want to calculate the index.
*/

periodFromTotalPeriod[ft, totalPeriod, index] :=
{
   N = length[ft]

   if index == 0  // Period is infinite.
      return undef

   if index <= N/2
      return totalPeriod / index
   else
      return -totalPeriod / (N - index)     
}


/** Calculate the period of the specified index given:
    
    [ft, samplePeriod, index]

    where:
      ft: the transformed sequence
      samplePeriod:  the period of a single sample (e.g. 10 ms)
      index:  the index for which we want to calculate the index.
*/

periodFromSamplePeriod[ft, samplePeriod, index] :=
{
   N = length[ft]

   totalPeriod = samplePeriod * N  // Is this valid with padded FFTs?

   println["totalPeriod is $totalPeriod, samplePeriod is $samplePeriod"]

   return periodFromTotalPeriod[ft, totalPeriod, index]
}

/** Calculate the frequency of the specified index given:
    
    [ft, samplePeriod, index]

    where:
      ft: the transformed sequence
      samplePeriod:  the period of each sample (e.g. 10 ms)
      index:  the index for which we want to calculate the index.

    Note that this gives the frequency.  If you want the *angular* frequency
    (omega) such as to pass to a sin or cos function, you need to multiply this
    frequency by 2 pi.
*/

frequencyFromSamplePeriod[ft, samplePeriod, index] :=
{
   if index == 0  // DC component
      return 0 * 1/samplePeriod  // Get the units right?

   println["length[ft] is " + length[ft]]

   return 1/periodFromSamplePeriod[ft, samplePeriod, index]
}

/** Plots a Fourier Transform as a magnitude-phase diagram. */
plotFourier[ft, showDC=true, centered=true, realOnly=false] :=
{
   N = length[ft]
   g = new graphics

   start = 0
   if showDC==false
      start = 1
   
   end = N-1
   if realOnly
      end = N div 2

   offset = 0
   modulus = N
   if centered
      offset = N div 2

   g.color[0, 0, 0, .7]
   for i = start to end
   {
      v = ft@i
      mag = magnitude[v]
      if realOnly
         arg = 2 Re[v]  // Real transform has complex conjugates at i and N-i
      else
         arg = arctan[Im[v], Re[v]]
      
      x = (i+offset) mod modulus
      g.color[0, 1, 0, 1]
      g.fillRectSides[x, -Re[arg], x+1, 0]
      g.color[0, 0, 0, .7]
      g.fillRectSides[x, -Re[mag], x+1, 0]
   }

   g.show[]
}


Download or view 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 19945 days, 13 hours, 48 minutes ago.