working on real optimization

fft.py

@@ -1,10 +1,12 @@
 #!/usr/local/bin/python2.3
 import math
+import sys
+
 pi=math.pi
 e=math.e
 j=complex(0,1)
 
-def fft(f):
+def fft(f,inv):
     n=len(f)
     if n==1:
         return f
@@ -22,7 +24,7 @@ def fft(f):
     for q in range(p): # 0,1
         fp = f[q::p]
         #print fp,'<< fp'
-        Fp = fft( fp )
+        Fp = fft( fp ,inv)
         Fout.extend( Fp )
 
     for u in range(m):
@@ -32,37 +34,78 @@ def fft(f):
             k = q1*m + u  # indices to Fout above that became scratch
             Fout[ k ] = scratch[0] # cuz e**0==1 in loop below
             for q in range(1,p):
-                t = e ** ( j*2*pi*k*q/n )
+                if inv:
+                    t = e ** ( j*2*pi*k*q/n )
+                else:                    
+                    t = e ** ( -j*2*pi*k*q/n )
                 Fout[ k ] += scratch[q] * t
 
     return Fout
 
-def real_fft( f ):
+def real_fft( f,inv ):
+    if inv:
+        sys.stderr.write('not impl\n')
+        sys.exit(1)
+
     N = len(f) / 2
 
     res = f[::2]
     ims = f[1::2]
 
     fp = [ complex(r,i) for r,i in zip(res,ims) ]
-
-    Fp = fft( fp )
+    Fp = fft( fp ,0 )
 
     F = []
     for k in range(N):
-        s2 = ( Fp[k] + Fp[-k] ) * .5
-        d2 = ( Fp[k] - Fp[-k] ).conjugate() * .5
-        F1k = complex( s2.real , d2.imag )
-        F2k = complex( s2.imag , d2.real )
-        F.append( F1k +  e ** ( -j*pi*k/N ) * F2k )
+        F1k = Fp[k] + Fp[-k].conjugate()
+        F2k = -j*(Fp[k] - Fp[-k].conjugate())
 
-    F.append( complex( Fp[0].real - Fp[0].imag , 0 ) ) 
+        F.append( ( F1k + e ** ( -j*pi*k/N ) * F2k ) * .5 )
+
+    F.append( complex( Fp[0].real - Fp[0].imag , 0 ) )
     return F
 
 def test(f=range(1024),ntimes=10):
     import time
     t0 = time.time()
     for i in range(ntimes):
-        fft(f)
+        fft(f,0)
     t1 = time.time()
     print '%ss'% (t1-t0)
 
+def main():
+    #fft_func = fft
+    fft_func = real_fft
+
+    tvec = [0.309655,0.815653,0.768570,0.591841,0.404767,0.637617,0.007803,0.012665]
+    Ftvec = [ complex(r,i) for r,i in zip(
+                [3.548571,-0.378761,-0.061950,0.188537,-0.566981,0.188537,-0.061950,-0.378761],
+                [0.000000,-1.296198,-0.848764,0.225337,0.000000,-0.225337,0.848764,1.296198] ) ]
+
+    F = fft_func( tvec,0 )
+
+    nerrs= 0
+    for i in range(len(Ftvec)/2 + 1):
+        if abs( F[i] - Ftvec[i] )> 1e-5:
+            print 'F[%d]: %s != %s' % (i,F[i],Ftvec[i])
+            nerrs += 1
+
+    print '%d errors in forward fft' % nerrs
+    if nerrs:
+        return
+
+    trec = fft_func( F , 1 )
+
+    for i in range(len(trec) ):
+        trec[i] /= len(trec)
+
+    for i in range(len(tvec) ):
+        if abs( trec[i] - tvec[i] )> 1e-5:
+            print 't[%d]: %s != %s' % (i,tvec[i],trec[i])
+            nerrs += 1
+
+    print '%d errors in reverse fft' % nerrs
+
+    
+if __name__ == "__main__":
+    main()