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
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)
3 loops, best of 3: 285 ms per loop
%timeit -r 3 -n 3 dft_py(ricker,dt,nf,df)
3 loops, best of 3: 1min 14s per loop
import numba
dft_nb=numba.autojit(dft_np)
fb=dft_nb(ricker,dt,nf,df)
plot(freq,abs(fb))
[<matplotlib.lines.Line2D at 0x11659ded0>]
%timeit -r 3 -n 3 dft_nb(ricker,dt,nf,df)
3 loops, best of 3: 1min 14s per loop
%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)
3 loops, best of 3: 2.82 s per loop
%%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
!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__
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)) .
fw=dft_fortran.dft_fv(ricker,dt,nf,df)
plot(freq,abs(fw))
[<matplotlib.lines.Line2D at 0x116b6b890>]
fw=dft_fortran.dft_fvec(ricker,dt,nf,df)
plot(freq,abs(fw))
[<matplotlib.lines.Line2D at 0x116bef210>]
%timeit -r 3 -n 3 dft_fortran.dft_fv(ricker,dt,nf,df)
3 loops, best of 3: 224 ms per loop
%timeit -r 3 -n 3 dft_fortran.dft_fvec(ricker,dt,nf,df)
3 loops, best of 3: 219 ms per loop