In [1]:

    

dt=0.001
nt=10000
tmax=nt*dt
time=np.linspace(0,tmax,nt)

df=1./tmax
fmax=100.
nf=int(fmax/df)
freq=np.linspace(0,fmax,nf)

print tmax,fmax


    

    




10.0 100.0

In [2]:

    

from scipy import signal
ricker=signal.ricker(nt,20.0)
plot(time,ricker,'k-')
xlabel("Time (s)",fontsize='large')
ylabel("Amplitude",fontsize='large')
savefig('ricker.eps')


    

In [3]:

    

def dft_py(w,dt,nf,df):
    fw=np.zeros(nf,dtype=np.complex64)
    t=np.linspace(0,w.shape[0]*dt,w.shape[0])
    for ifreq in xrange(nf):
        omega=2.*np.pi*ifreq*df
        fw[ifreq]=0.0j
        for it in xrange(w.shape[0]):
            fw[ifreq]+=w[it]*np.exp(-1.0j*omega*t[it])
        fw[ifreq]*=dt
    return fw

def dft_np(w,dt,nf,df):
    fw=np.zeros(nf,dtype=np.complex128)
    t=np.linspace(0,w.shape[0]*dt,w.shape[0])
    for ifreq in xrange(nf):
        omega=2.*np.pi*ifreq*df
        fw[ifreq]=np.dot(w,np.exp(-1.0j*omega*t))*dt
    return fw


    

In [4]:

    

fricker=dft_np(ricker,dt,nf,df)
plot(freq,abs(fricker),'k-')
xlabel('Frequency (Hz)',fontsize='large')
ylabel('Amplitude',fontsize='large')
savefig('spectrum.eps')


    

In [6]:

    

%timeit -r 3 -n 3 dft_np(ricker,dt,nf,df)


    

    




3 loops, best of 3: 285 ms per loop

In [7]:

    

%timeit -r 3 -n 3 dft_py(ricker,dt,nf,df)


    

    




3 loops, best of 3: 1min 14s per loop

In [11]:

    

import numba

dft_nb=numba.autojit(dft_np)


    

In [12]:

    

fb=dft_nb(ricker,dt,nf,df)
plot(freq,abs(fb))


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0x11659ded0>]






    

In [14]:

    

%timeit -r 3 -n 3 dft_nb(ricker,dt,nf,df)


    

    




3 loops, best of 3: 1min 14s per loop

In [1]:

    

%load_ext cythonmagic


    

In [2]:

    

%%cython
import cython
import numpy as np
cimport numpy as np
import cmath

ctypedef np.double_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def dft_cy(np.ndarray[double] w,double dt,int nf,double df):
    cdef np.ndarray[np.complex128_t] fw=np.zeros(nf,dtype=np.complex128)
    cdef np.ndarray[double,ndim=1] t=np.linspace(0.,w.shape[0]*dt,w.shape[0])
    cdef int ifreq,it
    cdef double pi=np.pi
    cdef DTYPE_t omega
    for ifreq in xrange(nf):
        omega=2.*pi*ifreq*df
        fw[ifreq]=0.0j
        for it in xrange(w.shape[0]):
            fw[ifreq]=fw[ifreq]+w[it]*cmath.exp(-1.0j*omega*t[it])
        fw[ifreq]=fw[ifreq]*dt
    return fw


    

In [18]:

    

%timeit -r 3 -n 3 dft_cy(ricker,dt,nf,df)


    

    




3 loops, best of 3: 2.82 s per loop

In [20]:

    

%%file dft_f.f90

subroutine dft_fv(fw,w,nt,dt,nf,df)
integer,intent(in):: nt,nf
real(kind=8),intent(in):: dt,df
real(kind=8),intent(in):: w(nt)
complex(kind=8),intent(out):: fw(nf)
integer it,ifreq
real(kind=8),parameter:: pi=3.141592654d0
real(kind=8):: omega,t(nt)
do it=1,nt
    t(it)=(it-1)*dt
enddo
do ifreq=1,nf
    omega=2.d0*pi*(ifreq-1)*df
    fw(ifreq)=dcmplx(0.d0,0.d0)
    do it=1,nt
        fw(ifreq)=fw(ifreq)+w(it)*cdexp(-dcmplx(0.d0,1.d0)*omega*t(it))
    enddo
    fw(ifreq)=fw(ifreq)*dt
enddo
end subroutine

subroutine dft_fvec(fw,w,nt,dt,nf,df)
integer,intent(in):: nt,nf
real(kind=8),intent(in):: dt,df
real(kind=8),intent(in):: w(nt)
complex(kind=8),intent(out):: fw(nf)
integer it,ifreq
real(kind=8),parameter:: pi=3.141592654d0
real(kind=8):: omega,t(nt)
do it=1,nt
    t(it)=(it-1)*dt
enddo
do ifreq=1,nf
    omega=2.d0*pi*(ifreq-1)*df
    fw(ifreq)=dot_product(w(:),cdexp(-dcmplx(0.d0,1.d0)*omega*t(:)))*dt
enddo
end subroutine


    

    




Overwriting dft_f.f90

In [21]:

    

!f2py -c -m dft_fortran dft_f.f90 --f90exec=/opt/local/bin/gfortran-mp-4.9 > log.txt


    

In [22]:

    

import dft_fortran
print dft_fortran.__doc__


    

    




This module 'dft_fortran' is auto-generated with f2py (version:2).
Functions:
  fw = dft_fv(w,dt,nf,df,nt=len(w))
  fw = dft_fvec(w,dt,nf,df,nt=len(w))
.

In [23]:

    

fw=dft_fortran.dft_fv(ricker,dt,nf,df)
plot(freq,abs(fw))


    

    
Out[23]:





[<matplotlib.lines.Line2D at 0x116b6b890>]






    

In [24]:

    

fw=dft_fortran.dft_fvec(ricker,dt,nf,df)
plot(freq,abs(fw))


    

    
Out[24]:





[<matplotlib.lines.Line2D at 0x116bef210>]






    

In [25]:

    

%timeit -r 3 -n 3 dft_fortran.dft_fv(ricker,dt,nf,df)


    

    




3 loops, best of 3: 224 ms per loop

In [26]:

    

%timeit -r 3 -n 3 dft_fortran.dft_fvec(ricker,dt,nf,df)


    

    




3 loops, best of 3: 219 ms per loop

In [26]: