# plot_fft.py
import numpy
from numpy.fft import fft
from numpy.fft import ifft
from numpy.matlib import real
from numpy import array
from numpy import zeros
import pylab
import math

print "plot_fft.py running"

nfft = 8; # signal amplitude, Inphase, no Quadrature
a=[0.0,0.7071,1.0,0.7071,0.0,-0.7071,-1.0,-0.7071]
print "a="
print a

b=fft(a)
print "b=fft(a)  complex  ="
print b
bb=abs(b)
for i in range(nfft): # normalize
  bb[i] = bb[i]/float(nfft)
  
print bb

rr = zeros(nfft/2)
xr = zeros(nfft/2)
rr[0]=0.0
xr[0]=0.0
for i in range(1,nfft/2):
  ar= b[i].real+b[nfft-i].real
  br= b[i].imag-b[nfft-i].imag
  rr[i] = math.sqrt(ar*ar+br*br)/float(nfft)
  xr[i] = float(i)

print 'rr=', rr

c=ifft(b)
print "c=inverse fft ifft(b)  complex ="
print c

d=real(c)
print "d=real(c) just real part"
print d

e=abs(c)
for i in range(nfft):
  e[i] = e[i]/float(nfft)
  
print "e=abs(c) normalized absolute value"
print e

x=array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0])
print "x="
print  x

print "generating plot_ffts_py.png"
pylab.plot(x, a)
pylab.xlabel("sample number, time")
pylab.ylabel("wave amplitude")
pylab.grid(True)
pylab.title("8 sample wave amplitude")
pylab.savefig("plot_ffts_py") # writes .png
pylab.show()

print "generating plot_fft_py.png"
pylab.plot(x, bb)
pylab.xlabel("frequency bin")
pylab.ylabel("abs spectrum")
pylab.grid(True)
pylab.title("FFT spectrum")
pylab.savefig("plot_fft_py") # writes .png
pylab.show()

pylab.plot(xr, rr)
pylab.xlabel("frequency bin")
pylab.ylabel("abs folded spectrum")
pylab.grid(True)
pylab.title("FFT unambiguous spectrum")
pylab.savefig("plot_fftu_py") # writes .png
pylab.show()

nfft = 16;
a16=[0.0, 0.3826, 0.7071, 0.9238, 1.0, 0.9238, 0.7071, 0.3826,
     0.0,-0.3826,-0.7071,-0.9238,-1.0,-0.9238,-0.7071,-0.3826]
print "a16="
print a16
x16=array([float(i) for i in range(nfft)])

b16=fft(a16)
print "b16=fft(a16)  complex  ="
print b16
bb16=abs(b16)
for i in range(nfft):
  bb16[i] = bb16[i]/float(nfft)
  
print "bb16 normalized absolute value of spectrum"
print bb16

rr16 = zeros(nfft/2)
xr16 = zeros(nfft/2)
rr16[0]=0.0
xr16[0]=0.0
for i in range(1,nfft/2):
  ar= b16[i].real+b16[nfft-i].real
  br= b16[i].imag-b16[nfft-i].imag
  rr16[i] = math.sqrt(ar*ar+br*br)/float(nfft)
  xr16[i] = float(i)


print "generating plot_fft16s_py.png"
pylab.plot(x16, a16)
pylab.xlabel("time sample")
pylab.ylabel("wave amplitude")
pylab.grid(True)
pylab.title("16 sample wave amplitude")
pylab.savefig("plot_fft16s_py") # writes .png
pylab.show()

print "generating plot_fft16_py.png"
pylab.plot(x16, bb16)
pylab.xlabel("frequency bin")
pylab.ylabel("abs spectrum")
pylab.grid(True)
pylab.title("FFT16 spectrum")
pylab.savefig("plot_fft16_py") # writes .png
pylab.show()

pylab.plot(xr16, rr16)
pylab.xlabel("frequency bin")
pylab.ylabel("abs folded spectrum")
pylab.grid(True)
pylab.title("FFT16 unambib spectrum")
pylab.savefig("plot_fft16u_py") # writes .png
pylab.show()

nfft = 128
f = open('tdoa_cwfm1.dat','r')
print f
aa=zeros(nfft)
xx=zeros(nfft)
for i in range(nfft):
  line = f.readline()
  print 'line=', line
  xx[i] = i
  aa[i] = float(line)
  print 'aa[',i,']=', aa[i]

f.close()

pylab.plot(xx, aa)
pylab.xlabel("time")
pylab.ylabel("amplitude")
pylab.title("tdoa_cwfm1.dat")
pylab.savefig("tdoa_cwfm1_dat.png") # writes .png
pylab.show()


cc=fft(aa)
cabs=abs(cc)
for i in range(nfft):
  cabs[i] = cabs[i]/float(nfft)
  
print cabs

pylab.plot(xx, cabs)
pylab.title("tdoa_cwfm1.dat spectrum")
pylab.xlabel("frequency bin")
pylab.ylabel("abs spectrum")
pylab.savefig("tdoa_cwfm1_spec.png") # writes .png
pylab.show()

rrr = zeros(nfft/2)
xrr = zeros(nfft/2)
rrr[0]=0.0
xrr[0]=0.0
for i in range(1,nfft/2):
  arr= cc[i].real+cc[nfft-i].real
  brr= cc[i].imag-cc[nfft-i].imag
  rrr[i] = math.sqrt(arr*arr+brr*brr)/float(nfft)
  xrr[i] = float(i)

pylab.plot(xrr, rrr)
pylab.title("tdoa_cwfm1.dat abs unambig spectrum")
pylab.xlabel("frequency bin")
pylab.ylabel("abs spectrum")
pylab.savefig("tdoa_cwfm1_specu.png") # writes .png
pylab.show()