fourierTemp.frink

//use englewoodarray.frink

a = new array
for line = lines["file:hourly.txt"]
a.push[eval[line]]
println[length[a]]

// Transform time-domain data into frequency-domain data
f = DFT[a]

g = new graphics

// fprime will contain just selected frequency components.
fprime = new array[]

len = length[f]
zeroUnit = 0 f@0
threshold = .1
vscale = 20

c = 0
str = ""

for y = f
{
x = (c + len/2) mod len      // Draw low-freq components together in middle

// Draw phases
g.color[0,1,0]
phase = arctan[Im[y],Re[y]]
g.fillRectSides[x,0,x+1,-5 phase]

// Draw magnitudes
g.color[0,0,0,.7]
mag = abs[y]
if (c != 0)                  // Don't draw DC component because it's huge.
g.fillRectSides[x,0,x+1,-vscale mag]

if (mag > threshold)                 // Decide which components to keep.
{

daysphase = (phase/(2 pi) * len)
mag2 = 2 mag
println["Used component \$c, mag is \$mag, phase is \$daysphase days"]
if c == 0
str = format[mag2/2, 1, 6]
if c > 0 && c < len/2
str = str + " +\n   " + format[mag2,1,6] + " cos[2 pi \$c/\$len d - " + format[phase,1,6] + "]"
fprime.push[y]
} else
fprime.push[zeroUnit]

c = c + 1
}

// Draw threshold
g.color[0,0,1,.5]
g.line[0,-vscale threshold,len,-vscale threshold]

println[str]
g.show[]
g.write["fourierEnglewood.svg",800,600]

// Turn the frequency-domain data back into a time-domain waveform
t = InverseDFT[fprime]

// Plot the time-domain data that was generated from the truncated
// frequency domain.
g2 = new graphics
c = 0
scale = 3
for y = t
{
g2.fillRectCenter[c, -scale Re[y], 1, scale]
c = c+1
}

// Plot original data
g2.color[0,0,1,.7]
c = 0
for y = a
{
g2.fillRectCenter[c, -scale Re[y], 1, scale]
c = c+1
}

// Test that sin approximation gives the same values
g2.color[1,0,0,.5]
for d = 0 to len-1
g2.fillRectCenter[d, -scale eval[str], 1, scale]

g2.show[]
g2.write["recreatedEnglewood.svg", 800, 600]

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 18321 days, 8 hours, 26 minutes ago.