Download or view FourierUtils.frink in plain text format
/** This file contains utilities for working with Fourier transforms, including
plotting the transform and figuring out what frequency components a
specific element of the Fourier transform represents. */
/** 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
dc = abs[ft@0]
if dc == 0
dc = 1
for i = start to end
{
v = ft@i
mag = abs[v]
if realOnly
arg = 0 // Real transform has complex conjugates at i and N-i
else
arg = arctan[Im[v], Re[v]]
x = (i+offset) mod modulus
transparency = clamp[mag/dc, 0 , 1]^(1/4)
// println["mag = $mag, dc = $dc, transparency = $transparency"]
g.color[0, 1, 0, transparency]
g.fillRectSides[x-1/2, -arg, x+1/2, 0]
g.color[0, 0, 0, .7]
g.fillRectSides[x-1/2, -Re[mag], x+1/2, 0]
}
return g
}
/** Adds a sinewave to the specified array. If the array is undef, this
creates a new array and returns it.
TODO: Specify phase?
*/
addSinewaveRaw[f, amplitude, startIndex, endIndex, array=undef] :=
{
if array == undef
array = new array[[endIndex+1],0]
len = endIndex-startIndex
for k=0 to len
{
prev = array@(k+startIndex)
if prev == undef
prev = 0
array@(k+startIndex) = prev + amplitude sin[2 pi f k]
}
return array
}
Download or view FourierUtils.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 19966 days, 13 hours, 33 minutes ago.