untitled
Random data to fit with a fourier transform
data fit with the FFT
1 2 3
4 5 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 29 var('n,x')
30 def wfunc(t):
31 return e^(i*((2*pi)*t))
32 33 def basisfunc(x,y,z):
34 return wfunc((x*y)/z)
35
36 37 def performFFT(spacetime_domain):
38 var('n,x')
39 print "spacetime_domain data"
40 print spacetime_domain
41 n = len(spacetime_domain)
42 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 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 58 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 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 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 93 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