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

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

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

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

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

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

import numba

dft_nb=numba.autojit(dft_np)

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

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

%load_ext cythonmagic

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

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

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

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

import dft_fortran
print dft_fortran.__doc__

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

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

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

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