Random data to fit with a fourier transform

data fit with the FFT

    1 # some Discrete and Fast fourier transforms in sage following
    2 # a Gilbert Strang Mit Linear Algebra lecture.
    3 
    4 #http://www.youtube.com/watch?v=M0Sa8fLOajA
    5 #not sure why strang's and everybody else's versions are negative/conjugate of each other, must be a change in convention...
    6 
    7 spacetime_domain = [100*random() for x in range(10)]
    8 
    9 def makeDFTNorm(n, w):
   10     rows = []
   11     for i in range(n):
   12         row = []
   13         for j in range(n):
   14             row.append(w^(-i*j))
   15         rows.append(row)
   16     return (1/sqrt(n))*matrix(rows)    
   17 
   18 def makeDFT(n, w):
   19     rows = []
   20     for i in range(n):
   21         row = []
   22         for j in range(n):
   23             row.append(w^(-i*j))
   24         rows.append(row)
   25     return matrix(rows)    
   26 
   27        
   28 #basis vectors  
   29 var('n,x')
   30 def wfunc(t):
   31     return e^(i*((2*pi)*t))
   32 #fbasis(x, y, z)=(w((-1)*(1/(x*y))*z));#*(1/z)
   33 def basisfunc(x,y,z):
   34     return wfunc((x*y)/z)
   35 
   36 # sage fast dft
   37 def performFFT(spacetime_domain):
   38     var('n,x')
   39     print "spacetime_domain data"
   40     print spacetime_domain
   41     n = len(spacetime_domain)
   42     #frequency_domain=vector(map(conjugate, vector(CDF, spacetime_domain).fft().list())).transpose();
   43     frequency_domain=vector(CDF, spacetime_domain).fft().transpose()
   44     print "real frequency domain "
   45     print map(lambda x: round(x.real_part(), 10), frequency_domain.list())
   46     #invdft = vector(map(conjugate, vector(map(conjugate, frequency_domain.list())).fft(direction='backward'))).transpose()
   47     invdft = vector(frequency_domain.list()).fft(direction='backward').transpose()
   48     reconstructed = map(real_part, invdft.list())
   49     plota = points([ (k, spacetime_domain[k]) for k in range(len(spacetime_domain)) ], rgbcolor=(1, 0, 0))
   50     plotb = points([ (k, reconstructed[k]) for k in range(len(reconstructed)) ], rgbcolor=(0, 1, 1))
   51     print "perfect point reconstruction"
   52     show(plota + plotb)
   53     var('axis')
   54     print axis
   55     basisVector = vector([basisfunc(axis,y,n) for y in range(n)])
   56     print basisVector.transpose()  
   57     #simplify/round frequency_domain for speedy plotting
   58     #frequency_domain = vector([round(frequency_domain.list()[l].real_part(), 10) + i*round(frequency_domain.list()[l].imag_part(), 10) for l in range(len(frequency_domain.list()))]).transpose()
   59     allBasesWithWeights=(basisVector*frequency_domain*(1/n)).list()[0].real_part()
   60     allBasesNoWeights=sum(basisVector*(1/n)).real_part()
   61     plotc = plot(allBasesWithWeights, xmin=0, xmax=n, color='red')
   62     plotd = plot(allBasesNoWeights, xmin=0, xmax=n, color='green')
   63     print "perfect point curve fitting"
   64     show(plota+plotb+plotc+plotd)
   65     var('n,x')
   66 
   67 
   68 performFFT(spacetime_domain)
   69 
   70 # strang slow DFT
   71 def performDFT(spacetime_domain):
   72     var('n,x')
   73     print "spacetime_domain data"
   74     print spacetime_domain
   75     n=len(spacetime_domain)
   76     myDFT=makeDFT(n, wfunc(1/n))
   77     frequency_domain=myDFT*(vector(CDF, spacetime_domain).transpose())
   78     print "real frequency domain "
   79     print map(lambda x: round(x.real_part(), 10), frequency_domain.list())
   80     myDFTH = Matrix(CDF, n,n, [x.conjugate() for x in myDFT.transpose().list()])
   81     invdft=(myDFTH*frequency_domain*(1/n))
   82     reconstructed = map(lambda x: round(x.real_part(), 10), invdft.list())
   83     print "perfect point reconstruction"
   84     #print reconstructed
   85     plota = points([ (k, spacetime_domain[k]) for k in range(len(spacetime_domain)) ], rgbcolor=(1, 0, 0))
   86     plotb = points([ (k, reconstructed[k]) for k in range(len(reconstructed)) ], rgbcolor=(0, 1, 1))
   87     show(plota + plotb)
   88     var('axis')
   89     print axis
   90     basisVector = vector([basisfunc(axis,y,n) for y in range(n)])
   91     print basisVector.transpose()
   92     #simplify/round frequency_domain for speedy plotting
   93     #frequency_domain = vector([(round(frequency_domain.list()[x].real_part(), 10) + i*round(frequency_domain.list()[x].imag_part(), 10)) for x in range(len(frequency_domain.list()))]).transpose()    
   94     allBasesWithWeights=(basisVector*frequency_domain*(1/n)).list()[0].real_part()
   95     allBasesNoWeights=sum(basisVector*(1/n)).real_part()
   96     plotc = plot(allBasesWithWeights, xmin=0, xmax=n, color='red')
   97     plotd = plot(allBasesNoWeights, xmin=0, xmax=n, color='green')
   98     print "perfect point curve fitting"
   99     show(plota+plotb+plotc+plotd)
  100     var('n,x')
  101 
  102 #performDFT(spacetime_domain)