2D Fourier Transform

with an extra gpu array for the 33th complex values



In [5]:

    
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 [6]:

    
#To put images inside the notebook

%matplotlib inline

Loading FFT routines



In [7]:

    
gridDIM = 64

size = gridDIM*gridDIM

axes0 = 0
axes1 = 1

makeC2C = 0
makeR2C = 1
makeC2R = 1

axesSplit_0 = 0
axesSplit_1 = 1

m = size

segment_axes0 = 0
segment_axes1 = 0

DIR_BASE = "/home/robert/Documents/new1/FFT/code/"

# FAFT 
_faft128_2D = ctypes.cdll.LoadLibrary( DIR_BASE+'FAFT128_2D_R2C.so' )
_faft128_2D.FAFT128_2D_R2C.restype = int
_faft128_2D.FAFT128_2D_R2C.argtypes = [ctypes.c_void_p, ctypes.c_void_p, 
                                         ctypes.c_float, ctypes.c_float, ctypes.c_int, 
                                         ctypes.c_int, ctypes.c_int, ctypes.c_int]

cuda_faft = _faft128_2D.FAFT128_2D_R2C

# Inv FAFT
_ifaft128_2D = ctypes.cdll.LoadLibrary( DIR_BASE+'IFAFT128_2D_C2R.so' )
_ifaft128_2D.IFAFT128_2D_C2R.restype = int
_ifaft128_2D.IFAFT128_2D_C2R.argtypes = [ctypes.c_void_p, ctypes.c_void_p, 
                                           ctypes.c_float, ctypes.c_float, ctypes.c_int, 
                                           ctypes.c_int, ctypes.c_int, ctypes.c_int]

cuda_ifaft = _ifaft128_2D.IFAFT128_2D_C2R

def fftGaussian(p,sigma): return np.exp( - p**2*sigma**2/2. )

Initializing Data

Gaussian



In [8]:

    
def Gaussian(x,sigma):
    return np.exp( - x**2/sigma**2/2.  )/(sigma*np.sqrt( 2*np.pi  ))

def fftGaussian(p,sigma):
    return np.exp( - p**2*sigma**2/2.  )



In [9]:

    
# Gaussian parameters
mu = 0
sigma = 1.

# Grid parameters
x_amplitude = 5.
p_amplitude = 6.                # With the traditional method p amplitude is fixed to: 2 * np.pi /( 2*x_amplitude ) 

dx = 2*x_amplitude/float(gridDIM)  # This is dx in Bailey's paper
dp = 2*p_amplitude/float(gridDIM)  # This is gamma in Bailey's paper

delta = dx*dp/(2*np.pi)

x_range = np.linspace( -x_amplitude, x_amplitude-dx, gridDIM) 
p = np.linspace( -p_amplitude, p_amplitude-dp, gridDIM) 

x = x_range[ np.newaxis, : ] 
y = x_range[ :, np.newaxis ] 

f =  Gaussian(x,sigma)*Gaussian(y,sigma)

plt.imshow( f, extent=[-x_amplitude , x_amplitude-dx, -x_amplitude , x_amplitude-dx] )

axis_font = {'size':'24'}
plt.text( 0., 5.1, '$W$' , **axis_font)
plt.colorbar()

#plt.ylim(0,0.44)


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 '    min = ', np.min(f)









    



 Amplitude x =  5.0
 Amplitude p =  6.0
        
sigma =  1.0
n     =  64
dx    =  0.15625
dp    =  0.1875
           standard fft dp =  0.628318530718      
    
delta =  0.0046627424734
    
The Gaussian extends to the numerical error in single precision:
    min =  2.21033491549e-12

$W$ TRANSFORM FROM AXES-0

After the transfom, f_gpu[:32, :] contains real values and f_gpu[32:, :] contains imaginary values. g33_gpu contains the 33th. complex values



In [10]:

    
f33 = np.zeros( [1 ,64], dtype = np.complex64 )



In [11]:

    
# One gpu array.

f_gpu = gpuarray.to_gpu( np.ascontiguousarray( f , dtype = np.float32 ) )
f33_gpu = gpuarray.to_gpu( np.ascontiguousarray( f33 , dtype = np.complex64 ) )

Forward Transform



In [12]:

    
# Executing FFT

cuda_faft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, delta, segment_axes0, axes0, makeR2C, axesSplit_0 )
cuda_faft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, delta, segment_axes1, axes1, makeC2C, axesSplit_0 )









    Out[12]:





1



In [13]:

    
plt.imshow( np.append( f_gpu.get()[:32, :], f33_gpu.get().real, axis=0 )/float(np.sqrt(size)), 
           extent=[-p_amplitude , p_amplitude-dp, 0, p_amplitude-dp] )

plt.colorbar()

axis_font = {'size':'24'}
plt.text( 0., 6.2, '$Re \\mathcal{F}(W)$',  **axis_font )

plt.xlim(-p_amplitude , p_amplitude-dp)
plt.ylim(0 , p_amplitude)









    Out[13]:





(0, 6.0)



In [14]:

    
plt.imshow( np.append( f_gpu.get()[32:, :], f33_gpu.get().imag, axis=0 )/float(np.sqrt(size)), 
           extent=[-p_amplitude , p_amplitude-dp, 0, p_amplitude-dp] )

plt.colorbar()

axis_font = {'size':'24'}
plt.text( 0., 6.2, '$Im \\mathcal{F}(W)$',  **axis_font )

plt.xlim(-p_amplitude , p_amplitude-dp)
plt.ylim(0 , p_amplitude)









    Out[14]:





(0, 6.0)

Central Section $p_x =0$



In [15]:

    
plt.figure(figsize=(10,10))

plt.plot( p[:33], np.append( f_gpu.get()[:32, 32], f33_gpu.get()[:, 32].real, axis=0 )/(f_gpu.get().size) , 
         'o', label='numerical')

plt.plot( p[:33], 4*fftGaussian(p,sigma)[:33] , label = 'analytical')

plt.legend(loc='upper left')

#plt.ylim(0,1.1)

plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[15]:





<matplotlib.text.Text at 0x7f5fd32f2990>



In [16]:

    
plt.figure(figsize=(10,10))

plt.semilogy( p[:33], np.append( f_gpu.get()[:32, 32], f33_gpu.get()[:, 32].real, axis=0 )/ (f_gpu.get().size), 
             'o', label='numerical')

plt.semilogy( p[:33], 4*fftGaussian(p,sigma)[:33] , label = 'analytical')

plt.legend(loc='upper left')

plt.ylim(0,1.1)

plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[16]:





<matplotlib.text.Text at 0x7f5fd5c0e950>



In [17]:

    
plt.plot( p[:33], np.append(f_gpu.get()[32:, 32], f33_gpu.get()[:, 32].imag, axis=0)/(f_gpu.get().size), 'o',
         label='numerical imag error')

plt.legend(loc='upper left')

#plt.ylim(0,1.1)

plt.xlabel('$p$',**axis_font)









    Out[17]:





<matplotlib.text.Text at 0x7f5fd339dd50>

Inverse Transform



In [18]:

    
# Executing iFFT

cuda_ifaft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, -delta, segment_axes1, axes1, makeC2C, axesSplit_0 )
cuda_ifaft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, -delta, segment_axes0, axes0, makeC2R, axesSplit_0 )









    Out[18]:





1



In [19]:

    
plt.imshow( f_gpu.get()/(float(size*size))  ,
           extent=[-x_amplitude , x_amplitude-dx, -x_amplitude, x_amplitude-dx] )

plt.colorbar()

axis_font = {'size':'24'}
plt.text( -1, 6.2, '$W$',  **axis_font )

plt.xlim(-x_amplitude , x_amplitude-dx)
plt.ylim(-x_amplitude , x_amplitude-dx)









    Out[19]:





(-5.0, 4.84375)



In [20]:

    
plt.figure(figsize=(10,10))

plt.plot( y, f_gpu.get()[:, 32]/(f_gpu.get().size) /(43.8*f_gpu.get().size), 
         'o', label='numerical')

plt.plot( y, 4*Gaussian(y, sigma) , label = 'analytical')

plt.legend(loc='upper left')

#plt.ylim(0,1.1)

plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[20]:





<matplotlib.text.Text at 0x7f5fd2a1add0>



In [21]:

    
# LOG PLOT 

plt.figure(figsize=(10,10))

plt.semilogy( y, f_gpu.get()[:, 32]/(f_gpu.get().size) /(450*f_gpu.get().size)  , 'o', label='numerical')

plt.semilogy( y, f[:, 32]  , '.', label='original')

plt.semilogy( y, Gaussian(y, sigma)/2.5 , label = 'analytical')

plt.legend(loc='upper left')

plt.ylim(1e-7,0.5)

plt.ylabel('$e^{- \\frac{x^2}{2 \\sigma^2 } }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[21]:





<matplotlib.text.Text at 0x7f5fd2b832d0>

$W$ TRANSFORM FROM AXES-1

After the transfom, f_gpu[:, :32] contains real values and f_gpu[:, 32:] contains imaginary values. f33_gpu contains the 33th. complex values



In [22]:

    
f33 = np.zeros( [64, 1], dtype = np.complex64 )



In [23]:

    
# One gpu array.

f_gpu = gpuarray.to_gpu( np.ascontiguousarray( f , dtype = np.float32 ) )
f33_gpu = gpuarray.to_gpu( np.ascontiguousarray( f33 , dtype = np.complex64 ) )

Forward Transform



In [24]:

    
# Executing FFT

cuda_faft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, delta, segment_axes1, axes1, makeR2C, axesSplit_1 )
cuda_faft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, delta, segment_axes0, axes0, makeC2C, axesSplit_1 )









    Out[24]:





1



In [25]:

    
plt.imshow( np.append( f_gpu.get()[:, :32], f33_gpu.get().real, axis=1 )/float(np.sqrt(size)), 
           extent=[-p_amplitude , 0, -p_amplitude, p_amplitude-dp] )

plt.colorbar()

axis_font = {'size':'24'}
plt.text( -3.0, 6.2, '$Re \\mathcal{F}(W)$',  **axis_font )

plt.xlim(-p_amplitude , 0)
plt.ylim(-p_amplitude , p_amplitude-dp)









    Out[25]:





(-6.0, 5.8125)



In [26]:

    
plt.imshow( np.append( f_gpu.get()[:, 32:], f33_gpu.get().imag, axis=1 )/float(np.sqrt(size)), 
           extent=[-p_amplitude , 0, -p_amplitude, p_amplitude-dp] )

plt.colorbar()

axis_font = {'size':'24'}
plt.text( -3.0, 6.2, '$Im \\mathcal{F}(W)$',  **axis_font )

plt.xlim(-p_amplitude , 0)
plt.ylim(-p_amplitude , p_amplitude-dp)









    Out[26]:





(-6.0, 5.8125)

Central Section $p_y =0$



In [27]:

    
plt.figure(figsize=(10,10))

plt.plot( p[:33], np.append( f_gpu.get()[32, :32], f33_gpu.get()[32, :].real, axis=0 )/(f_gpu.get().size) , 
         'o', label='numerical')

plt.plot( p[:33], 4*fftGaussian(p,sigma)[:33] , label = 'analytical')

plt.legend(loc='upper left')

#plt.ylim(0,1.1)

plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[27]:





<matplotlib.text.Text at 0x7f5fd25b5110>



In [28]:

    
plt.figure(figsize=(10,10))

plt.semilogy( p[:33], np.append( f_gpu.get()[32, :32], f33_gpu.get()[32, :].real, axis=0 )/ (f_gpu.get().size), 
             'o', label='numerical')

plt.semilogy( p[:33], 4*fftGaussian(p,sigma)[:33] , label = 'analytical')

plt.legend(loc='upper left')

plt.ylim(0,1.1)

plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[28]:





<matplotlib.text.Text at 0x7f5fd284ae90>



In [29]:

    
plt.plot( p[:33], np.append(f_gpu.get()[32, 32:], f33_gpu.get()[32, :].imag, axis=0)/(f_gpu.get().size), 'o',
         label='numerical imag error')

plt.legend(loc='upper left')

#plt.ylim(0,1.1)

plt.xlabel('$p$',**axis_font)









    Out[29]:





<matplotlib.text.Text at 0x7f5fd25e5990>

Inverse Transform



In [30]:

    
# Executing iFFT

cuda_ifaft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, -delta, segment_axes0, axes0, makeC2C, axesSplit_1 )
cuda_ifaft( int(f_gpu.gpudata), int(f33_gpu.gpudata), dx, -delta, segment_axes1, axes1, makeC2R, axesSplit_1 )









    Out[30]:





1



In [31]:

    
plt.imshow( f_gpu.get()/float(size)  ,
           extent=[-x_amplitude , x_amplitude-dx, -x_amplitude , x_amplitude-dx] )

plt.colorbar()

axis_font = {'size':'24'}
plt.text( -1, 5.2, '$W$',  **axis_font )

plt.xlim(-x_amplitude , x_amplitude-dx)
plt.ylim(-x_amplitude , x_amplitude-dx)









    Out[31]:





(-5.0, 4.84375)



In [32]:

    
plt.figure(figsize=(10,10))

plt.plot( y, f_gpu.get()[32, :]/(f_gpu.get().size) /(43.8*f_gpu.get().size), 
         'o', label='numerical')

plt.plot( y, 4*Gaussian(y, sigma) , label = 'analytical')

plt.legend(loc='upper left')

#plt.ylim(0,1.1)

plt.ylabel('$e^{- \\frac{\\sigma x^2}{2} }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[32]:





<matplotlib.text.Text at 0x7f5fd23a07d0>



In [33]:

    
# LOG PLOT 

plt.figure(figsize=(10,10))

plt.semilogy( y, f_gpu.get()[32, :]/(f_gpu.get().size) /(450*f_gpu.get().size)  , 'o', label='numerical')

plt.semilogy( y, f[32, :]  , '.', label='original')

plt.semilogy( y, Gaussian(y, sigma)/2.5 , label = 'analytical')

plt.legend(loc='upper left')

plt.ylim(1e-7,0.5)

plt.ylabel('$e^{- \\frac{x^2}{2 \\sigma^2 } }$',**axis_font)
plt.xlabel('$p$',**axis_font)









    Out[33]:





<matplotlib.text.Text at 0x7f5fd22edb90>



In [ ]: