In [85]:

    

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

    

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


    

    
Out[86]:





<matplotlib.text.Text at 0x110fc5290>






    

In [87]:

    

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

    

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


    

    
Out[88]:





<matplotlib.text.Text at 0x1113d9e10>






    

In [90]:

    

ncut=100
%timeit -r 3 -n 3 dft_np(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 9.38 ms per loop

In [91]:

    

ncut=1000
%timeit -r 3 -n 3 dft_np(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 33.7 ms per loop

In [92]:

    

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


    

    




3 loops, best of 3: 285 ms per loop

In [93]:

    

ncut=100
%timeit -r 3 -n 3 dft_py(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 784 ms per loop

In [94]:

    

ncut=1000
%timeit -r 3 -n 3 dft_py(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 7.84 s per loop

In [95]:

    

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


    

    




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

In [96]:

    

import numba

dft_nb=numba.autojit(dft_py)


    

In [99]:

    

ncut=100
%timeit -r 3 -n 3 dft_nb(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 786 ms per loop

In [100]:

    

ncut=1000
%timeit -r 3 -n 3 dft_nb(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 7.87 s per loop

In [101]:

    

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


    

    




3 loops, best of 3: 1min 18s 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
ctypedef np.complex128_t DTYPEC_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 [105]:

    

ncut=100
%timeit -r 3 -n 3 dft_cy(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 29.3 ms per loop

In [106]:

    

ncut=1000
%timeit -r 3 -n 3 dft_cy(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 284 ms per loop

In [107]:

    

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


    

    




3 loops, best of 3: 2.86 s per loop

In [109]:

    

%%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 [110]:

    

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


    

In [111]:

    

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

    

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


    

    
Out[112]:





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






    

In [113]:

    

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


    

    
Out[113]:





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






    

In [114]:

    

ncut=100
%timeit -r 3 -n 3 dft_fortran.dft_fv(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 2.09 ms per loop

In [115]:

    

ncut=1000
%timeit -r 3 -n 3 dft_fortran.dft_fv(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 21.4 ms per loop

In [116]:

    

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


    

    




3 loops, best of 3: 218 ms per loop

In [116]:

    




    

In [117]:

    

ncut=100
%timeit -r 3 -n 3 dft_fortran.dft_fvec(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 2.24 ms per loop

In [118]:

    

ncut=1000
%timeit -r 3 -n 3 dft_fortran.dft_fvec(ricker[:ncut],dt,nf,df)


    

    




3 loops, best of 3: 21.2 ms per loop

In [119]:

    

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


    

    




3 loops, best of 3: 216 ms per loop