In [1]:
import numpy as np
import ctypes
from ctypes import *
import pycuda.gpuarray as gpuarray
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import math
In [2]:
%matplotlib inline
In [3]:
gridDIM = 2048
size = gridDIM
axes0 = 0
m = size
DIR_BASE = "/home/robert/Documents/new1/FFT/code/FAFT4096/"
# FAFT
_faft4096_1D = ctypes.cdll.LoadLibrary( DIR_BASE+'FAFT4096_1D_R2C.so' )
_faft4096_1D.FAFT4096_1D_R2C.restype = int
_faft4096_1D.FAFT4096_1D_R2C.argtypes = [ctypes.c_void_p, ctypes.c_void_p,
ctypes.c_float, ctypes.c_float, ctypes.c_int]
cuda_faft = _faft4096_1D.FAFT4096_1D_R2C
# IFAFT
_ifaft4096_1D = ctypes.cdll.LoadLibrary( DIR_BASE+'IFAFT4096_1D_C2R.so' )
_ifaft4096_1D.IFAFT4096_1D_C2R.restype = int
_ifaft4096_1D.IFAFT4096_1D_C2R.argtypes = [ctypes.c_void_p, ctypes.c_void_p,
ctypes.c_float, ctypes.c_float, ctypes.c_int]
cuda_ifaft = _ifaft4096_1D.IFAFT4096_1D_C2R
In [4]:
def fftGaussian(p,sigma):
return np.exp( - p**2*sigma**2/2. )
In [5]:
def Gaussian(x,mu,sigma):
return np.exp( -(x-mu)**2/(2*sigma**2) )/np.sqrt( np.sqrt(np.pi)*sigma )
def fftGaussian(p,mu,sigma):
return np.exp( - p**2*sigma**2/2. )*np.exp( -1j*mu*p ) * np.sqrt( sigma/np.sqrt(np.pi) )
In [20]:
# Gaussian parameters
mu = 1.5
#mu = 0
sigma = 1.
# Grid parameters
x_amplitude = 8.
p_amplitude = 5. # With the traditional method p amplitude is fixed to: 2 * np.pi /( 2*x_amplitude )
dx = 2*x_amplitude/float(size) # This is beta in Bailey's paper
dp = 2*p_amplitude/float(size) # This is gamma in Bailey's paper
delta = dx*dp/(2*np.pi)
x = np.linspace( -x_amplitude, x_amplitude-dx, size)
p = np.linspace( -p_amplitude, p_amplitude-dx, size)
f = Gaussian(x, mu, sigma)
plt.plot(x, f,'-')
axis_font = {'size':'24'}
plt.text( -5, 1.1, 'Full Signal on the Host (2048 points)' , **axis_font)
plt.ylabel('$e^{-\\frac{1}{\\sigma }x^2}$',**axis_font)
plt.xlabel('$x$',**axis_font)
plt.ylim(0,1)
print ' Amplitude x = ',x_amplitude
print ' Amplitude p = ',p_amplitude
print ' '
print 'sigma = ', sigma
print 'n = ', x.size
print 'dx = ', dx
print 'dp = ', dp
print ' standard fft dp = ',2 * np.pi /( 2*x_amplitude ) , ' '
print ' '
print 'delta = ', delta
print ' '
print 'The Gaussian extends to the numerical error in single precision:'
print ' extremes ', f[0], ' , ', f[-1]
In [7]:
f1025 = np.zeros( [ 1 ], dtype = np.complex64 )
In [8]:
# Copy data to GPU
segment = 0
f_gpu = gpuarray.to_gpu( np.ascontiguousarray( f , dtype = np.float32 ) )
f1025_gpu = gpuarray.to_gpu( np.ascontiguousarray( f1025 , dtype = np.complex64 ) )
In [9]:
# LOG PLOT
plt.figure(figsize=(10,10))
plt.semilogy( x, f_gpu.get() , '-', label='numerical', markersize=1)
plt.semilogy( x, Gaussian(x, mu, sigma) , 'r--' ,label = 'analytical')
plt.legend(loc='upper left')
plt.ylim(1e-8,1)
plt.ylabel('$e^{- \\frac{x^2}{2 \\sigma^2 } }$',**axis_font)
plt.xlabel('$x$',**axis_font)
Out[9]:
In [10]:
# Executing FFT
cuda_faft( int(f_gpu.gpudata), int(f1025_gpu.gpudata), dx, delta, segment )
Out[10]:
In [11]:
# Normalization
norm = np.sum(np.abs(f_gpu.get())**2) + np.sum(np.abs(f1025_gpu.get())**2)
norm *= dp
norm = np.sqrt(norm)
f_gpu /= norm
f1025_gpu /= norm
In [12]:
plt.figure(figsize=(10,10))
plt.plot( f_gpu.get() , '.', label='numerical', markersize=1)
plt.legend(loc='upper left')
#plt.ylim(-0.3 , 1.1)
plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)
Out[12]:
In [13]:
def ReconstructFFT(f):
n = f.shape[0]
freal_half = f_gpu.get()[:n/2]
freal = np.append( freal_half , f1025_gpu.get().real )
freal = np.append( freal , freal_half[:0:-1] )
fimag_half = f_gpu.get()[n/2:]
fimag = np.append( fimag_half , np.array([0.]) )
fimag = np.append( fimag , -fimag_half[:0:-1] )
return freal + 1j*fimag
In [14]:
plt.figure(figsize=(10,10))
plt.plot( ReconstructFFT( f_gpu.get() ).real , '-' )
plt.plot( ReconstructFFT( f_gpu.get() ).real , '--' , label='Real')
plt.plot( ReconstructFFT( f_gpu.get() ).imag , '-' )
plt.plot( ReconstructFFT( f_gpu.get() ).imag , '--' ,label='Imag')
plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)
plt.legend(loc='upper left')
Out[14]:
In [15]:
# Standard fft :
# Observe how much memory is lost in zeros
plt.figure(figsize=(10,10))
plt.plot( np.fft.fftshift( np.fft.fft( np.fft.fftshift(f) ) ).real ,'-', label='Real')
plt.plot( np.fft.fftshift( np.fft.fft( np.fft.fftshift(f) ) ).real ,'--' )
plt.plot( np.fft.fftshift( np.fft.fft( np.fft.fftshift(f) ) ).imag ,'-', label='Imag')
plt.plot( np.fft.fftshift( np.fft.fft( np.fft.fftshift(f) ) ).imag ,'--')
plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)
plt.legend(loc='upper left')
Out[15]:
In [16]:
# For Inverse Transform, enter "-delta"
cuda_ifaft( int(f_gpu.gpudata), int(f1025_gpu.gpudata), dx, -delta, segment )
Out[16]:
In [17]:
# Normalization
norm = np.sum(np.abs(f_gpu.get())**2) + np.sum(np.abs(f1025_gpu.get())**2)
norm *= dx
norm = np.sqrt(norm)
f_gpu /= norm
f1025_gpu /= norm
In [18]:
plt.figure(figsize=(10,10))
plt.plot( x, f_gpu.get(), '-', label='numerical')
plt.plot( x, Gaussian(x, mu, sigma) , '--',label = 'analytical')
plt.legend(loc='upper left')
#plt.ylim(0,0.5)
plt.ylabel('$e^{- \\frac{x^2}{2 \\sigma^2 } }$',**axis_font)
plt.xlabel('$p$',**axis_font)
Out[18]:
In [19]:
# LOG PLOT
plt.figure(figsize=(10,10))
plt.semilogy( x, f_gpu.get() , '-', label='numerical')
plt.semilogy( x, f , '--', label='original')
plt.semilogy( x, Gaussian(x, mu, sigma) , label = 'analytical')
plt.legend(loc='upper left')
plt.ylim(1e-8,0.5)
plt.ylabel('$e^{- \\frac{x^2}{2 \\sigma^2 } }$',**axis_font)
plt.xlabel('$p$',**axis_font)
Out[19]:
In [19]: