# Fourier.frink

``` // 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[] } ```

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 18864 days, 22 hours, 54 minutes ago.