1D Fast Accurate Fourier Transform - 2048 points

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

segment = 0

DIR_BASE = "/home/rcabrera/Documents/source/Python/FFT-dev/FAFT2048_C2C-tmp6/"

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

cuda_faft = _faft4096_1D.FAFT4096_1D_C2C


    

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

    

# 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( -10, 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]


    

    




 Amplitude x =  8.0
 Amplitude p =  5.0
        
sigma =  1.0
n     =  2048
dx    =  0.0078125
dp    =  0.0048828125
           standard fft dp =  0.392699081699      
    
delta =  6.07127926223e-06
    
The Gaussian extends to the numerical error in single precision:
  extremes  1.89746614536e-20  ,  5.28788942035e-10






    

In [7]:

    

# Copy data to GPU

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


    

In [8]:

    

# 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)


    

    




/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:482: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)






    
Out[8]:





<matplotlib.text.Text at 0x7f13854837d0>






    

In [9]:

    

# Executing FFT

cuda_faft( int(f_gpu.gpudata), dx, delta, segment )


    

    
Out[9]:





1

In [10]:

    

# Normalization 
norm = np.sum(np.abs(f_gpu.get())**2)
norm *= dp
norm = np.sqrt(norm)

f_gpu    /= norm


    

In [11]:

    

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

plt.plot( f_gpu.get().real , '.', label='Real', markersize=1 )
plt.plot( f_gpu.get().imag , 'r-', label='Imag', 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[11]:





<matplotlib.text.Text at 0x7f1384f32790>






    

In [12]:

    

# 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[12]:





<matplotlib.legend.Legend at 0x7f1384e50cd0>






    

In [13]:

    

# For Inverse Transform, enter "-delta"

cuda_faft( int(f_gpu.gpudata), dx, -delta, segment )


    

    
Out[13]:





1

In [14]:

    

# Normalization 
norm = np.sum(np.abs(f_gpu.get())**2)
norm *= dx
norm = np.sqrt(norm)

f_gpu    /= norm


    

In [15]:

    

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

plt.plot( x, f_gpu.get().real, '-', 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[15]:





<matplotlib.text.Text at 0x7f1384f154d0>






    

In [16]:

    

# LOG PLOT 

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

plt.semilogy( x, f_gpu.get().real  , '-', 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,1)

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


    

    
Out[16]:





<matplotlib.text.Text at 0x7f138492ded0>






    

In [16]: