In [3]:

    

# menu column at left and images list at right. to be refined..
# plus some styling
from IPython.display import HTML
import urllib2
TOC="https://gist.githubusercontent.com/astyonax/1d7f49edefd9728b7da0/raw"
TOC=urllib2.urlopen(TOC).read()
INAV="https://gist.github.com/astyonax/815eacf44e2efd5eb0a7/raw"
INAV=urllib2.urlopen(INAV).read()
style="https://gist.githubusercontent.com/astyonax/d53a8eb1c556168a65f5/raw"
style=urllib2.urlopen(style).read()
HTML("""
<style>{2}</style>
<style>
    .container{{
        width: 90% !important;
        max-width:90% !important;
    }}
    #img_nav{{
    
        max-width:10% !important;
        padding: 1px !important;
        z-index:10000000000 !important;
    }}
    #thumb_nav{{
        max-width:100%;
    }}
    div.output_area img {{
    max-width:95% !important;
    }}
</style>

<div id="toc"></div>
<script type="text/javascript">
{0}
</script>
<div id="img_nav"></div>
<script type="text/javascript">
{1}
</script>""".format(TOC,INAV,style))


    

    
Out[3]:

In [1]:

    

%pylab inline
import numpy as np
import pylab as pl
import scipy
import numexpr as ne
from numba import jit
import sympy


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

def morlet(eta,w0):
    return ne.evaluate("3.141592653589793**-.25*exp(-.5 * eta**2)*exp(1j*w0*eta)")

def C(s,w0):
    return np.pi**-.25/np.sqrt(2*s)*exp(.25*(w0-sqrt(2+w0**2))**2)

def wlet(y,s,tau,dt,w0):
    
    out=y[0]*morlet((0*dt-tau)/s,w0)
    for i in xrange(1,y.shape[0]): # this is slow but, now, I do not know how to go faster
        out+=y[i]*morlet((i*dt-tau)/s,w0)
    return dt/np.sqrt(s)*out
    

def wlet_2(y,s,tau,dt,w0):
    
    ## THIS IS SUPER MEMORY HUNGRY!!!!
    NT=np.arange(y.shape[0])
    out=morlet((dt*NT[:,np.newaxis,np.newaxis]-tau)/s,w0)
    
    out=(y[:,np.newaxis,np.newaxis]*out).sum(axis=0)
    return dt/np.sqrt(s)*out



f2s=lambda f,w0:(w0+sqrt(2+w0**2))/(4*np.pi*f)
s2f=lambda s,w0:(w0+sqrt(2+w0**2))/(4*np.pi*s)

def wlet_trans(data,fmin,fmax,df,dt_sample,dt,w0=50):
    """
    This function performs the 1D wavelet transform of 'data'
    using the morlet wavelet.
    The code is based on the paper: http://dx.doi.org/10.1098/rsif.2014.0672
    
    Input:
    -data is a 1D numpy.ndarray
    -fmin,  fmax, and df define the frequency range spanned by the transform 
     as in np.arange(fmin,fmax+df,df)
    -dt_sample is the sampling rate of the input signal
    -dt is the transform time stepping
    -w0 -- see ref.. determines time vs spatial resolution
    
    Output:
    -spectrogram
    -times
    -frequencies
    """
    smax=f2s(fmin,w0)
    smin=f2s(fmax,w0)
    
    taus=np.arange(0,data.shape[0]*dt_sample,dt)
    ff=np.arange(fmin,fmax+df,df)
    ss=f2s(ff,w0)
    
    tauss,ss=np.meshgrid(taus,ss)
    X,Y=tauss.shape
    print "transform size:",X*Y,"bins"
    WL=wlet
    if X*Y<41000:
        WL=wlet_2
    transformed_data=wlet(data,s=ss,tau=tauss,dt=dt_sample,w0=w0)/C(ss,w0=w0)[::-1]
    #return C(ss,w0=w0),taus,ff
    print transformed_data.shape
    return transformed_data,taus,ff


    

In [17]:

    

etas=np.arange(-4,4,.001)
w0s=[5,10,20]
N=len(w0s)
f,ax=pl.subplots(1,N,figsize=(N*5,5),sharex='all',sharey='all')
for w0,ax in zip(w0s,ax):
    ax.plot(etas,morlet(etas,w0).real,lw=1)
    ax.plot(etas,morlet(etas,w0).imag,'--',lw=1)
    ax.set_title(r"$\omega_0={:d}$".format(w0))
    ax.set_xlabel(r"$\eta$")


    

In [18]:

    

# Test with synthetic signal
dt=1/500. # 500Hz
t=np.arange(0,3.5,dt)
data=np.sin(2*np.pi*30*t)*(t<1)
data+=np.sin(2*np.pi*25*t)*((t>1)&(t<2))
data+=np.sin(2*np.pi*t*20)*((t>2) & (t<3))
data+=np.sin(2*np.pi*t*22)*((t>2) & (t<3))
data+=np.sin(2*np.pi*t*40)*.3
pl.figure(figsize=(15,4))
pl.plot(t,data)

w0s=[30,40,50]
N=len(w0s)
f,ax=pl.subplots(1,N,figsize=(N*5,5),sharex='all',sharey='all')
for w0,ax in zip(w0s,ax):
    tdata,taus,ff=wlet_trans(data,15,60,1,dt,dt*2,w0)
    extend=[0,taus.max(),ff.min(),ff.max()]
    ax.imshow(np.abs(tdata)**.5,aspect="auto",cmap=pl.cm.gist_heat,extent=extend,interpolation="none",origin="bottom")
    ax.set_title(r"$\omega_0={:d}$".format(w0))
    ax.set_xlabel(r"$time$")


    

    




transform size: 40250 bins
(46, 875)
transform size: 40250 bins
(46, 875)
transform size: 40250 bins
(46, 875)






    













    

In [61]:

    

T=np.arange(0,800)
S=np.arange(1,30)
TT,SS=np.meshgrid(T,S)
print morlet((TT[:,np.newaxis]*dt-TT)/SS,10.).shape


    

    




(29, 29, 800)

In [2]:

    

import sympy
sympy.init_printing()


    

In [3]:

    

s,t,tau,f,w0,eta,y,psi=sympy.symbols("s t tau f w0 eta y psi",positive=1,real=1)

psi=(sympy.pi**(-sympy.Rational(1,4))*sympy.exp(-sympy.Rational(1,2) *eta**2)*sympy.exp(sympy.I*w0*eta))
psi=psi.subs(eta,(t-tau)/s)
y=sympy.sin(2*sympy.pi*f*t)
psi,y


    

    
Out[3]:





$$\left ( \frac{1}{\sqrt[4]{\pi}} e^{- \frac{\left(t - \tau\right)^{2}}{2 s^{2}}} e^{\frac{i w_{0}}{s} \left(t - \tau\right)}, \quad \sin{\left (2 \pi f t \right )}\right )$$

In [4]:

    

WL=sympy.integrate(psi*y,(t,-sympy.oo,sympy.oo))


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-4-ce86fdd258c3> in <module>()
----> 1 WL=sympy.integrate(psi*y,(t,-sympy.oo,sympy.oo))

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/utilities/decorator.pyc in threaded_func(expr, *args, **kwargs)
     33                                       func(expr.rhs, *args, **kwargs))
     34             else:
---> 35                 return func(expr, *args, **kwargs)
     36 
     37     return threaded_func

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/integrals.pyc in integrate(*args, **kwargs)
   1230     if isinstance(integral, Integral):
   1231         return integral.doit(deep=False, meijerg=meijerg, conds=conds,
-> 1232                              risch=risch, manual=manual)
   1233     else:
   1234         return integral

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/integrals.pyc in doit(self, **hints)
    465                 and not function.is_Poly and \
    466                     (xab[1].has(oo, -oo) or xab[2].has(oo, -oo)):
--> 467                 ret = try_meijerg(function, xab)
    468                 if ret is not None:
    469                     function = ret

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/integrals.pyc in try_meijerg(function, xab)
    442                     x, a, b = xab
    443                     try:
--> 444                         res = meijerint_definite(function, x, a, b)
    445                     except NotImplementedError:
    446                         from sympy.integrals.meijerint import _debug

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in meijerint_definite(f, x, a, b)
   1744                 _debug('  Non-real splitting point.')
   1745                 continue
-> 1746             res1 = _meijerint_definite_2(f.subs(x, x + c), x)
   1747             if res1 is None:
   1748                 _debug('  But could not compute first integral.')

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in _meijerint_definite_2(f, x)
   1853     for g, explanation in _guess_expansion(f, x):
   1854         _debug('Trying', explanation)
-> 1855         res = _meijerint_definite_3(g, x)
   1856         if res is not None and res[1] != False:
   1857             return res

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in _meijerint_definite_3(f, x)
   1870     if f.is_Add:
   1871         _debug('Expanding and evaluating all terms.')
-> 1872         ress = [_meijerint_definite_4(g, x) for g in f.args]
   1873         if all(r is not None for r in ress):
   1874             conds = []

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in _meijerint_definite_4(f, x, only_double)
   1921 
   1922     # Try two G functions.
-> 1923     gs = _rewrite2(f, x)
   1924     if gs is not None:
   1925         for full_pb in [False, True]:

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in _rewrite2(f, x)
   1573     for recursive in [False, True]:
   1574         for fac1, fac2 in l:
-> 1575             g1 = _rewrite_single(fac1, x, recursive)
   1576             g2 = _rewrite_single(fac2, x, recursive)
   1577             if g1 and g2:

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in _rewrite_single(f, x, recursive)
   1501     try:
   1502         F, strip, _ = mellin_transform(f, x, s, integrator=my_integrator,
-> 1503                                        simplify=False, needeval=True)
   1504         g = my_imt(F, s, x, strip)
   1505     except IntegralTransformError:

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/transforms.pyc in mellin_transform(f, x, s, **hints)
    351     hankel_transform, inverse_hankel_transform
    352     """
--> 353     return MellinTransform(f, x, s).doit(**hints)
    354 
    355 

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/transforms.pyc in doit(self, **hints)
    116             try:
    117                 return self._compute_transform(self.function,
--> 118                     self.function_variable, self.transform_variable, **hints)
    119             except IntegralTransformError:
    120                 pass

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/transforms.pyc in _compute_transform(self, f, x, s, **hints)
    293 
    294     def _compute_transform(self, f, x, s, **hints):
--> 295         return _mellin_transform(f, x, s, **hints)
    296 
    297     def _as_integral(self, f, x, s):

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/transforms.pyc in wrapper(*args, **kwargs)
    194         def wrapper(*args, **kwargs):
    195             noconds = kwargs.pop('noconds', default)
--> 196             res = func(*args, **kwargs)
    197             if noconds:
    198                 return res[0]

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/transforms.pyc in _mellin_transform(f, x, s_, integrator, simplify)
    218     # convergence of the integral.
    219     s = _dummy('s', 'mellin-transform', f)
--> 220     F = integrator(x**(s - 1) * f, x)
    221 
    222     if not F.has(Integral):

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/integrals/meijerint.pyc in my_integrator(f, x)
   1495         if r is not None:
   1496             res, cond = r
-> 1497             res = _my_unpolarify(hyperexpand(res, rewrite='nonrepsmall'))
   1498             return Piecewise((res, cond),
   1499                              (Integral(f, (x, 0, oo)), True))

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/hyperexpand.pyc in hyperexpand(f, allow_hyper, rewrite)
   2480         if not r.has(nan, zoo, oo, -oo):
   2481             return r
-> 2482     return f.replace(hyper, do_replace).replace(meijerg, do_meijer)

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in replace(self, query, value, map, simultaneous, exact)
   1349             return expr
   1350 
-> 1351         rv = bottom_up(self, rec_replace, atoms=True)
   1352 
   1353         # restore original expressions for Dummy symbols

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/simplify.pyc in bottom_up(rv, F, atoms, nonbasic)
   4077         if rv.args:
   4078             args = tuple([bottom_up(a, F, atoms, nonbasic)
-> 4079                 for a in rv.args])
   4080             if args != rv.args:
   4081                 rv = rv.func(*args)

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/simplify.pyc in bottom_up(rv, F, atoms, nonbasic)
   4077         if rv.args:
   4078             args = tuple([bottom_up(a, F, atoms, nonbasic)
-> 4079                 for a in rv.args])
   4080             if args != rv.args:
   4081                 rv = rv.func(*args)

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/simplify.pyc in bottom_up(rv, F, atoms, nonbasic)
   4080             if args != rv.args:
   4081                 rv = rv.func(*args)
-> 4082             rv = F(rv)
   4083         elif atoms:
   4084             rv = F(rv)

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in rec_replace(expr)
   1334             result = _query(expr)
   1335             if result or result == {}:
-> 1336                 new = _value(expr, result)
   1337                 if new is not None and new != expr:
   1338                     mapping[expr] = new

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in <lambda>(expr, result)
   1278                 _value = lambda expr, result: value(*expr.args)
   1279             elif callable(value):
-> 1280                 _value = lambda expr, result: value(*expr.args)
   1281             else:
   1282                 raise TypeError(

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/hyperexpand.pyc in do_meijer(ap, bq, z)
   2477     def do_meijer(ap, bq, z):
   2478         r = _meijergexpand(G_Function(ap[0], ap[1], bq[0], bq[1]), z,
-> 2479                            allow_hyper, rewrite=rewrite)
   2480         if not r.has(nan, zoo, oo, -oo):
   2481             return r

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/hyperexpand.pyc in _meijergexpand(func, z0, allow_hyper, rewrite)
   2376 
   2377     slater1 = powdenest(slater1.subs(z, z0), polar=True)
-> 2378     slater2 = powdenest(slater2.subs(t, 1/z0), polar=True)
   2379     if not isinstance(cond2, bool):
   2380         cond2 = cond2.subs(t, 1/z)

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/simplify.pyc in powdenest(eq, force, polar)
   2440 
   2441     if polar:
-> 2442         eq, rep = polarify(eq)
   2443         return unpolarify(powdenest(unpolarify(eq, exponents_only=True)), rep)
   2444 

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/simplify/simplify.pyc in polarify(eq, subs, lift)
   2168         return eq
   2169     reps = dict([(s, Dummy(s.name, polar=True)) for s in eq.free_symbols])
-> 2170     eq = eq.subs(reps)
   2171     return eq, dict([(r, s) for s, r in reps.items()])
   2172 

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in subs(self, *args, **kwargs)
    890             rv = self
    891             for old, new in sequence:
--> 892                 rv = rv._subs(old, new, **kwargs)
    893                 if not isinstance(rv, Basic):
    894                     break

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in _subs(self, old, new, **hints)
   1004         rv = self._eval_subs(old, new)
   1005         if rv is None:
-> 1006             rv = fallback(self, old, new)
   1007         return rv
   1008 

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in fallback(self, old, new)
    977                     continue
    978                 arg = arg._subs(old, new, **hints)
--> 979                 if not _aresame(arg, args[i]):
    980                     hit = True
    981                     args[i] = arg

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in _aresame(a, b)
   1672     """
   1673     from .function import AppliedUndef, UndefinedFunction as UndefFunc
-> 1674     for i, j in zip_longest(preorder_traversal(a), preorder_traversal(b)):
   1675         if i != j or type(i) != type(j):
   1676             if ((isinstance(i, UndefFunc) and isinstance(j, UndefFunc)) or

/home/astyonax/.anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in __iter__(self)
   1824         return next(self._pt)
   1825 
-> 1826     def __iter__(self):
   1827         return self
   1828 

KeyboardInterrupt: 

In [11]:

    

result=sympy.simplify(WL)
result


    

    
Out[11]:





$$\frac{\sqrt{2} \sqrt{\pi} s}{2 e^{\frac{f^{2} s^{2}}{2} + f s w_{0} + \frac{w_{0}^{2}}{2}}} \left(e^{2 f s w_{0}} \sin{\left (f \tau \right )} + i e^{2 f s w_{0}} \cos{\left (f \tau \right )} + \sin{\left (f \tau \right )} - i \cos{\left (f \tau \right )}\right)$$

In [16]:

    

sympy.simplify(result*result.conjugate())


    

    
Out[16]:





$$\frac{\pi s^{2}}{2 e^{f^{2} s^{2} + 2 f s w_{0} + w_{0}^{2}}} \left(e^{4 f s w_{0}} + 4 e^{2 f s w_{0}} \sin^{2}{\left (f \tau \right )} - 2 e^{2 f s w_{0}} + 1\right)$$

In [ ]: