Commits

mike fc committed 493e54a

tests for fftpack_lite now in separate file

Comments (0)

Files changed (2)

lib_pypy/numpypy/fft/fftpack_lite.py

 	# a = [r, r1+i1, r2+i2, ...] ==> [r, r1, i1, r2, i2, ...]
 
 	n = (len(wrk) - 15) // 2
-	print "guessed n:", n
 	res = np.zeros(n, dtype='f8')
 	res[0] = np.real(a[0])
 	if n % 2 != 0:
 		res[1:n-1:2] = np.real(a[1:len(a)-1])
 		res[2:n  :2] = np.imag(a[1:len(a)-1])
 		res[-1] = np.real(a[-1])
-	print "INV: input n", n, res
 	res_ptr = ffi.cast("double *", res.__array_interface__['data'][0])
 	wrk_ptr = ffi.cast("double *", wrk.__array_interface__['data'][0])
 	fftpack_c.rfftb(n, res_ptr, wrk_ptr)
 	res2 = res[0:n2-1:2] + 1j * res[1:n2  :2]
 	return res2 / n
 
-
-#-------------------------------------------------------------------
-# Helpers
-#-------------------------------------------------------------------
-def fft1(x):
-    L = len(x)
-    phase = -2j*np.pi*(np.arange(L)/float(L))
-    phase = np.arange(L).reshape(-1, 1) * phase
-    return np.sum(x*np.exp(phase), axis=1)
-
-def assert_array_almost_equal(a1, a2):
-	assert len(a1)==len(a2)
-	assert (np.absolute(a1 - a2)).sum() < 1e-6
-
-if __name__ == '__main__':
-	# a   = np.array([85.0, 76.0, 42.0, 26.0, 51.0, 40.0, 79.0, 30.0], dtype='f8')
-	# a   = np.array([85.0, 76.0, 42.0, 26.0, 51.0], dtype='f8')
-	# print "-"*80
-	# n = len(a)
-	# wrk = rffti(n)
-	# res = rfftf(a, wrk)
-	# print res
-	# print "-"*80
-	# # tres = fft1(a)
-	# # assert_array_almost_equal(res, tres[:4])
-	# res = rfftb(res, wrk)
-	# print res
-	# print a
-	a = np.array([ 1.+2.j,  2.+3.j,  2.-2.j,  3.-1.j], dtype='complex128')
-	print "ORG:", a
-	print "-"*80
-	n = len(a)
-	wrk = cffti(n)
-	res = cfftf(a, wrk)
-	print "FIN:", res
-	print "-"*80
-	wrk = cffti(n)
-	res = cfftb(res, wrk)
-	print "BAK:", res
-	print "ORG:", a
-
-
-
-
-
-
-
-
-
-
-

lib_pypy/numpypy/fft/test_fftpack_lite.py

+#!/usr/bin/env python
+
+import numpy as np
+from fftpack_lite import *
+
+#-------------------------------------------------------------------
+# Helpers
+#-------------------------------------------------------------------
+def fft1(x):
+    L = len(x)
+    phase = -2j*np.pi*(np.arange(L)/float(L))
+    phase = np.arange(L).reshape(-1, 1) * phase
+    return np.sum(x*np.exp(phase), axis=1)
+
+def assert_array_almost_equal(a1, a2):
+	assert len(a1)==len(a2), "%i != %i" % (len(a1), len(a2))
+	assert (np.absolute(a1 - a2)).sum() < 1e-6
+
+# real ffts
+for n in [8, 11, 17]:
+    print "Real", n,
+    a = np.random.random(n)
+    wrk = rffti(n)
+    res = rfftf(a, wrk)
+    # slow version of fft
+    tres = fft1(a)
+    t = n//2 + 1
+    assert_array_almost_equal(res, tres[:t])
+
+    # inverse fft
+    res = rfftb(res, wrk)
+    assert_array_almost_equal(res, a)
+    print "OK"
+
+# complex ffts
+for n in [8, 11, 17]:
+    print "Complex", n, 
+    a = np.random.random(n) + 1j * np.random.random(n)
+    wrk = cffti(n)
+    res = cfftf(a, wrk)
+    # slow version of fft
+    tres = fft1(a)
+    assert_array_almost_equal(res, tres)
+
+    # inverse fft
+    res = cfftb(res, wrk)
+    assert_array_almost_equal(res, a)
+    print 'OK'
+
+
+
+
+
+