thesis/fourier/fourier.py

   1 #!/usr/bin/python
   2
   3 import sys
   4 import os
   5 import cmath
   6 import math
   7 import copy
   8
   9 import threading, Queue
  10 import time
  11
  12 class Thread(threading.Thread):
  13         def __init__(self,x, inverse=False, spawn_threads=1):
  14                 threading.Thread.__init__(self)
  15                 self.x = x
  16                 self.inverse = inverse
  17                 self.X = []
  18                 self.spawn_threads = spawn_threads
  19
  20         def run(self):
  21                 #print "Thread " + str(threading.current_thread()) + " runs\n"
  22                 self.X = fft(self.x, self.inverse, self.spawn_threads)
  23
  24
  25 def dft(x, inverse = False):
  26
  27         #print "Using Brute Force"
  28
  29         N = len(x)
  30         inv = -1 if not inverse else 1
  31
  32         X = [0] * N
  33         for k in range(0, N):
  34                 for n in range(0, N):
  35                         X[k] += x[n] * cmath.exp(inv * 2.0j * cmath.pi * float(k * n) / float(N))
  36                 if inverse:
  37                         X[k] /= float(N)
  38
  39         return X
  40
  41 def fft(x, inverse = False, spawn_threads=1):
  42
  43         N = len(x)
  44
  45         if N == 1:
  46                 return [x[0]]
  47
  48         inv = -1 if not inverse else 1
  49
  50         if N % 2 != 0:
  51                 return dft(x, inverse)
  52
  53         #print "Active threads " + str(threading.active_count()) + "; in thread " + str(threading.current_thread())
  54         if threading.active_count() < spawn_threads:
  55                 threads = [Thread(x[0::2]), Thread(x[1::2])]
  56
  57                 for t in threads:
  58                         t.start()
  59
  60                 for t in threads:
  61                         t.join()
  62
  63                 X = threads[0].X + threads[1].X
  64         else:
  65                 X = fft(x[0::2]) + fft(x[1::2])
  66
  67
  68         for k in range(0, N/2):
  69                 t = X[k]
  70                 X[k] = t + cmath.exp(inv * 2.0j * cmath.pi * float(k) / float(N)) * X[k + (N/2)]
  71                 X[k+(N/2)] = t - cmath.exp(inv * 2.0j * cmath.pi * float(k) / float(N)) * X[k + (N/2)]
  72
  73         if not inverse:
  74                 return X
  75         else:
  76                 return map(lambda j : j / 2.0, X)
  77
  78 def transpose(lists, defval=0.0):
  79         if not lists: return []
  80         return map(lambda *row : [elem or defval for elem in row], *lists)
  81
  82
  83 def ifft(x):
  84         return fft(x, True)
  85
  86 def fft2d(x, inverse = False):
  87         xt = transpose(x)
  88         temp = [[]] * len(xt)
  89         for k in range(0, len(xt)):
  90                 temp[k] = fft(xt[k], inverse)
  91
  92         temp = transpose(temp)
  93         X = [[]] * len(temp)
  94         for k in range(0, len(temp)):
  95                 X[k] = fft(temp[k], inverse)
  96         return X
  97
  98 def offset2d(x):
  99
 100         for i in range(0, len(x)/2):
 101                 for j in range(0, len(x[i])/2):
 102                         t = x[i][j]
 103                         x[i][j] = x[i+len(x)/2][j+len(x[i])/2]
 104                         x[i+len(x)/2][j+len(x[i])/2] = t
 105
 106         for i in range(len(x)/2, len(x)):
 107                 for j in range(0, len(x[i])/2):
 108                         t = x[i][j]
 109                         x[i][j] = x[i-len(x)/2][j+len(x[i])/2]
 110                         x[i-len(x)/2][j+len(x[i])/2] = t
 111
 112
 113         return x
 114
 115 def main(argv):
 116         return 0
 117
 118 if __name__ == "__main__":
 119         sys.exit(main(sys.argv))