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

Loading FFT routines


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

Initializing a Gaussian


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]


 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]:
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]:
<matplotlib.text.Text at 0x7f43d3e0aa50>

Forward Transform


In [10]:
# Executing FFT

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


Out[10]:
1

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]:
<matplotlib.text.Text at 0x7f43d360d4d0>

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]:
<matplotlib.legend.Legend at 0x7f43d3526d50>

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]:
<matplotlib.legend.Legend at 0x7f43d30825d0>

Inverse


In [16]:
# For Inverse Transform, enter "-delta"

cuda_ifaft( int(f_gpu.gpudata), int(f1025_gpu.gpudata), dx, -delta, segment )


Out[16]:
1

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]:
<matplotlib.text.Text at 0x7f43d3537d50>

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]:
<matplotlib.text.Text at 0x7f43d2fdf310>

In [19]: