DIY¶

0. Import python libraries¶

import numpy as np
import pylab
from pylab import detrend_mean
import math

1.Choose and implement the wavelet function¶

def wave_bases(mother,k,scale,param):
  """Computes the wavelet function as a function of Fourier frequency
   used for the CWT in Fourier space (Torrence and Compo, 1998)
   -- This def is called automatically by def wavelet --

   _____________________________________________________________________
   Inputs:
          mother - a string equal to 'Morlet'
          k      - a vectorm the Fourier frequecies
          scale  - a number, the wavelet scale
          param  - the nondimensional parameter for the wavelet function

   Outputs:
          daughter       - a vector, the wavelet function
          fourier_factor - the ratio os Fourier period to scale
          coi            - a number, the cone-of-influence size at the scale
          dofmin         - a number, degrees of freedom for each point in the wavelet power (Morlet = 2)

   Call function:
          daughter,fourier_factor,coi,dofmin = wave_bases(mother,k,scale,param)
   _____________________________________________________________________

• expnt Check it out:

  expnt = -pow(scale*k-k0,2)/2*(k>0)
  

+---------------------------------------+
|               scale[a1] = 32          |
|             +-----------------+       |
|             |-0               |       |
|             |                 |       |
|len(k) = 512 |                 |       |
|             |                 |       |
|             |               -0|       |
|             +-----------------+       |
+---------------------------------------+

• norm Check it out:

  norm = math.sqrt(scale*k[1])*(pow(math.pi,-0.25))*math.sqrt(len(k))
  

+---------------------------------------+
|          scale[a1] = 32               |
|       +-----------------+             |
|       |2.66        39.07|             |
|       +-----------------+             |
|                                       |
+---------------------------------------+

• daughter Check it out:

  daughter = []                                           # define daughter as a list
  for ex in expnt:                                        # for each value scale (equal to next pow of 2)
          daughter.append(norm*math.exp(ex))
  k = np.array(k)                                         # turn k to array
  daughter = np.array(daughter)                           # transform in array
  daughter = daughter*(k>0)                               # Heaviside step function
  

+---------------------------------------+
|               scale[a1] = 32          |
|             +-----------------+       |
|             |0.0              |       |
|             |                 |       |
|len(k) = 512 |                 |       |
|             |                 |       |
|             |              0.0|       |
|             +-----------------+       |
+---------------------------------------+

2. Find the next pow of 2¶

def nextpow2(i):
    n = 2
    while n < i: n = n * 2
    return n

3. Compute wavelet power spectrum¶

def wavelet(Y,dt,mother,param):#,pad,dj,s0,J1,mother,param):
      """Computes the wavelet continuous transform of the vector Y, by definition:

         W(a,b) = sum(f(t)*psi[a,b](t) dt)               a dilate/contract
         psi[a,b](t) = 1/sqrt(a) psi(t-b/a)              b displace

         Only Morlet wavelet (k0=6) is used
         The wavelet basis is normalized to have total energy = 1 at all scales
         _____________________________________________________________________
         Input:
              Y - time series
              dt - sampling rate
              mother - the mother wavelet function
              param - the mother wavelet parameter

         Output:
              ondaleta - wavelet bases at scale 10 dt
              wave - wavelet transform of Y
              period - the vector of "Fourier"periods ( in time units) that correspond to the scales
              scale - the vector of scale indices, given by S0*2(j*DJ), j =0 ...J1
              coi - cone of influence

         Call function:
         ondaleta, wave, period, scale, coi = wavelet(Y,dt,mother,param)
         _____________________________________________________________________

n1 = len(Y)                                                     # time series length
s0 = 2*dt                                                       # smallest scale of the wavelet
dj = 0.25                                                       # spacing between discrete scales

• J1 Check it out

J1= int(np.floor((np.log10(n1*dt/s0))/np.log10(2)/dj))          # J1+1 total os scales

• Call nextpow2 Check it out

if (pad ==1) :
        base2 = nextpow2(n1)                                    #call det nextpow2
n = base2

• k Check it out:

  # construct wavenumber array used in transform
  # simetric eqn 5
  k = np.arange(n/2)
  import math
  k_pos,k_neg=[],[]
  for i in range(0,n/2+1):
          k_pos.append(i*((2*math.pi)/(n*dt)))                    # frequencies as in eqn5
          k_neg = k_pos[::-1]                                     # inversion vector
          k_neg = [e * (-1) for e in k_neg]                       # negative part
          k_neg = k_neg[1:-1]                                     # delete the first value of k_neg = last value of k_pos
  k = np.concatenate((k_pos,k_neg), axis =1)                      # vector of symmetric
  

+---------------------------------------+
|                 1                     |
|             +------+                  |
|             | 0.0  |                  |
|             |      |                  |
|len(k) = 512 |      |                  |
|             |      |                  |
|             |-0.049|                  |
|             +------+                  |
+---------------------------------------+

• f Check it out:

  # compute fft of the padded time series
  f = np.fft.fft(x,n)
  

+---------------------------------------+
|                 1                     |
|             +------------+            |
|             |-9.33e-15+0j|            |
|             |            |            |
|len(k) = 512 |            |            |
|             |            |            |
|             |31.9 -5.55j |            |
|             +------------+            |
+---------------------------------------+