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 19979 days, 9 hours, 55 minutes ago.