Commits

mike fc committed c80c68c

a replacement fftpack_lite that uses cffi. works in isolation with numpy, but no idea on how to get this into numpypy

Comments (0)

Files changed (3)

lib_pypy/numpypy/fft/Makefile

+all: fftpack.o
+	gcc -shared -W1,-soname,libfftpack.so.1 -o libfftpack.so.1.0 fftpack.o
+
+fftpack.o: fftpack.c
+	gcc -O2 -Wall -fPIC -c fftpack.c
+
+clean:
+	rm *.o libfft*
+

lib_pypy/numpypy/fft/fftpack.py

 version of the FFTPACK routines.
 
 """
-__all__ = ['fft','ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn',
-           'irfftn', 'rfft2', 'irfft2', 'fft2', 'ifft2', 'fftn', 'ifftn']
+#__all__ = ['fft','ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn',
+#           'irfftn', 'rfft2', 'irfft2', 'fft2', 'ifft2', 'fftn', 'ifftn']
+# Not supporting any of the multidimensional ffts just yet - mikefc.
+__all__ = ['fft','ifft', 'rfft', 'irfft', 'hfft', 'ihfft']
 
-from numpy.core import asarray, zeros, swapaxes, shape, conjugate, \
+from numpy import asarray, zeros, swapaxes, shape, conjugate, \
      take
 import fftpack_lite as fftpack
 

lib_pypy/numpypy/fft/fftpack_lite.py

+#!/usr/bin/env python
+
+import numpy as np
+
+# This is a replacement for numpy's fftpack_litemodule.c
+# written by mikefc for numpypy. targetting fftpack with cffi.
+
+#-------------------------------------------------------------------
+# CFFI interface to fftpack.c
+#-------------------------------------------------------------------
+import cffi
+ffi = cffi.FFI()
+
+fftpack_h = """
+void cfftf(int N, double data[], const double wrk[]);
+void cfftb(int N, double data[], const double wrk[]);
+void cffti(int N, double wrk[]);
+void rfftf(int N, double data[], const double wrk[]);
+void rfftb(int N, double data[], const double wrk[]);
+void rffti(int N, double wrk[]);
+"""
+
+ffi.cdef(fftpack_h)
+fftpack_c = ffi.dlopen("./libfftpack.so.1.0")
+
+#-------------------------------------------------------------------
+# 1D FFT of real values
+#-------------------------------------------------------------------
+def rffti(n):
+	"""Init function for real-valued fft
+	@param n length of fft to be performed.
+	@return appropriately sized working array 
+	"""
+	wrk = np.zeros(2 * n + 15, dtype='f8')
+	wrk_ptr = ffi.cast("double *", wrk.__array_interface__['data'][0])
+	fftpack_c.rffti(n, wrk_ptr)
+	return wrk
+
+def rfftf(a, wrk):
+	"""Perform an fft on 'a' assuming that it only contains real values
+	@param a real-valued 1D array
+	@wrk a previously initialised working array
+	@return complex ndarray 
+	"""
+	a = a.copy().astype('f8') # C code works in-place. Avoid overwriting
+	n = len(a)
+	a_ptr   = ffi.cast("double *",   a.__array_interface__['data'][0])
+	wrk_ptr = ffi.cast("double *", wrk.__array_interface__['data'][0])
+	fftpack_c.rfftf(n, a_ptr, wrk_ptr)
+	# The result from fftpack_c.rfftf needs to be unscrambled to create
+	# the complex array. Each pair of doubles represents a complex
+	# number, except for the first value.
+	# i.e. a is [r, r1, i1, r2, i2, ....]
+	#    res = [r, r1+i1, r2+i2, ... r]
+	if n % 2 == 0:
+		# Even length ffts have a real term at 0 and n//2
+		# the terms 1..n//2 are the conj of the terms n//2+1..n
+		res = np.zeros(n//2+1, dtype='complex128')
+		res[0]    = a[0]
+		res[1:n//2] = a[1:n-1:2] + 1j * a[2:n:2]
+		res[n//2] = a[-1]
+	else:
+		# Odd length ffts have a real term at 0
+		# r, c1, c2 ... cn, conj(cn), ... conj(c2), conj(c1)
+		res = np.zeros((n+1)//2, dtype='complex128')
+		res[0]    = a[0]
+		res[1:(n+1)//2] = a[1:n-1:2] + 1j * a[2:n:2]
+	return res
+
+def rfftb(a, wrk):
+	# reconfigure input for fftpack
+	# 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  :2] = np.real(a[1:len(a)])
+		res[2:n  :2] = np.imag(a[1:len(a)])
+	else:
+		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)
+	return res / n
+
+#-------------------------------------------------------------------
+# 1D FFT of general (complex) values
+#-------------------------------------------------------------------
+def cffti(n):
+	"""Init function for real-valued fft
+	@param n length of fft to be performed.
+	@return appropriately sized working array 
+	"""
+	wrk = np.zeros(4 * n + 15, dtype='f8')
+	wrk_ptr = ffi.cast("double *", wrk.__array_interface__['data'][0])
+	fftpack_c.cffti(n, wrk_ptr)
+	return wrk
+
+def cfftf(a, wrk):
+	""" Compute the general forwards complex fft """
+	n = len(a)  # n complex numbers
+	n2 = 2 * n  # represented by 2*n doubles
+	
+	# unpack complex values for c
+	# [r1+i1, r2+i2, .... rn+in]
+	#   ==>  [r1, i1, r2, i2, .... rn, in]
+	res = np.zeros(n2, dtype='f8')
+	res[0:n2-1:2] = np.real(a)
+	res[1:n2  :2] = np.imag(a)
+
+	res_ptr = ffi.cast("double *", res.__array_interface__['data'][0])
+	wrk_ptr = ffi.cast("double *", wrk.__array_interface__['data'][0])
+	fftpack_c.cfftf(n, res_ptr, wrk_ptr)
+
+	# pack complex values for python
+	# [r1, i1, r2, i2, .... rn, in]
+	#   ==>  [r1+i1, r2+i2, .... rn+in]
+	res2 = np.zeros(n, dtype='complex128')
+	res2 = res[0:n2-1:2] + 1j * res[1:n2  :2]
+	return res2
+
+def cfftb(a, wrk):
+	"""Compute the general backwards complex fft
+	Identical to cfftf except for fftpack_c function call
+	"""
+	n = len(a)  # n complex numbers
+	n2 = 2 * n  # represented by 2*n doubles
+
+	# unpack complex values for c
+	# [r1+i1, r2+i2, .... rn+in]
+	#   ==>  [r1, i1, r2, i2, .... rn, in]
+	res = np.zeros(n2, dtype='f8')
+	res[0:n2-1:2] = np.real(a)
+	res[1:n2  :2] = np.imag(a)
+
+	res_ptr = ffi.cast("double *", res.__array_interface__['data'][0])
+	wrk_ptr = ffi.cast("double *", wrk.__array_interface__['data'][0])
+	fftpack_c.cfftb(n, res_ptr, wrk_ptr)
+
+	# pack complex values for python
+	# [r1, i1, r2, i2, .... rn, in]
+	#   ==>  [r1+i1, r2+i2, .... rn+in]
+	res2 = np.zeros(n, dtype='complex128')
+	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
+
+
+
+
+
+
+
+
+
+
+