//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]