Commit before breaking everything

[matches/honours.git] / research / analysis and stuff / fourier / fourier.py

diff --git a/research/analysis and stuff/fourier/fourier.py b/research/analysis and stuff/fourier/fourier.py
new file mode 100644 (file)
index 0000000..893a3f8
--- /dev/null
+++ b/research/analysis and stuff/fourier/fourier.py
@@ -0,0 +1,119 @@
+#!/usr/bin/python
+
+import sys
+import os
+import cmath
+import math
+import copy
+
+import threading, Queue
+import time
+
+class Thread(threading.Thread):
+       def __init__(self,x, inverse=False, spawn_threads=1):
+               threading.Thread.__init__(self)
+               self.x = x
+               self.inverse = inverse
+               self.X = []
+               self.spawn_threads = spawn_threads
+               
+       def run(self):
+               #print "Thread " + str(threading.current_thread()) + " runs\n"
+               self.X = fft(self.x, self.inverse, self.spawn_threads)
+               
+
+def dft(x, inverse = False):
+
+       #print "Using Brute Force"
+
+       N = len(x)
+       inv = -1 if not inverse else 1
+       
+       X = [0] * N
+       for k in range(0, N):
+               for n in range(0, N):
+                       X[k] += x[n] * cmath.exp(inv * 2.0j * cmath.pi * float(k * n) / float(N))
+               if inverse:
+                       X[k] /= float(N)
+       
+       return X
+
+def fft(x, inverse = False, spawn_threads=1):
+
+       N = len(x)
+
+       if N == 1:
+               return [x[0]]
+
+       inv = -1 if not inverse else 1
+       
+       if N % 2 != 0:
+               return dft(x, inverse)
+
+       #print "Active threads " + str(threading.active_count()) + "; in thread " + str(threading.current_thread())
+       if threading.active_count() < spawn_threads:
+               threads = [Thread(x[0::2]), Thread(x[1::2])]
+       
+               for t in threads:
+                       t.start()
+
+               for t in threads:
+                       t.join()
+               
+               X = threads[0].X + threads[1].X
+       else:
+               X = fft(x[0::2]) + fft(x[1::2])
+       
+       
+       for k in range(0, N/2):
+               t = X[k]
+               X[k] = t + cmath.exp(inv * 2.0j * cmath.pi * float(k) / float(N)) * X[k + (N/2)]
+               X[k+(N/2)] = t - cmath.exp(inv * 2.0j * cmath.pi * float(k) / float(N)) * X[k + (N/2)]
+
+       if not inverse:
+               return X
+       else:
+               return map(lambda j : j / 2.0, X)
+
+def transpose(lists, defval=0.0):
+       if not lists: return []
+       return map(lambda *row : [elem or defval for elem in row], *lists)
+
+
+def ifft(x):
+       return fft(x, True)
+
+def fft2d(x, inverse = False):
+       xt = transpose(x)
+       temp = [[]] * len(xt)
+       for k in range(0, len(xt)):
+               temp[k] = fft(xt[k], inverse)
+       
+       temp = transpose(temp)
+       X = [[]] * len(temp)
+       for k in range(0, len(temp)):
+               X[k] = fft(temp[k], inverse)
+       return X
+
+def offset2d(x):
+
+       for i in range(0, len(x)/2):
+               for j in range(0, len(x[i])/2):
+                       t = x[i][j]
+                       x[i][j] = x[i+len(x)/2][j+len(x[i])/2]
+                       x[i+len(x)/2][j+len(x[i])/2] = t 
+
+       for i in range(len(x)/2, len(x)):
+               for j in range(0, len(x[i])/2):
+                       t = x[i][j]
+                       x[i][j] = x[i-len(x)/2][j+len(x[i])/2]
+                       x[i-len(x)/2][j+len(x[i])/2] = t 
+
+
+       return x
+
+def main(argv):
+       return 0
+
+if __name__ == "__main__":
+       sys.exit(main(sys.argv))