Example for sine waves¶

First written on Nov,19, 2014:
Author: tmiyama at gmail

Introduction¶

Wavelet anaysis tool by Torrence and Compo [1998] is very popular in climatology. However, Liu et al. [2007] has pointed out that the method by Torrence and Compo [1998] has a bias in favor of large scales.

This notebook is the translation of Liu et al.'s matlab test program
http://ocgweb.marine.usf.edu/~liu/wavelet_test_ElNino3_YLiu.m
to python.

This notebook use an artificial sine curve for a test correponding to Liu et al. [2007].

Notebook for NINO3 SST is here: http://nbviewer.ipython.org/urls/dl.dropbox.com/s/ryeoz3wxrkg1hhz/wavelet_test_ElNino3_Liu.ipynb?dl=0

References:

Here is the comment in the original matlab script.

% Rectification of the bias in the Wavelet power spectrum with the data set
% (Nino3.dat) given by Torrence and Compo (1998). This code is modified
% from wavetest.m, the example script provided by Torrence and Compo % (1998), to demonstrate how to rectify the bias in wavelet power spectrum.
% This code generates Figure 4 of Liu et al. (2007).
%
% Yonggang Liu, 2006.4.12
%
% E-mail: [email protected]
% http://ocgweb.marine.usf.edu/~liu/wavelet.html
% % References:
% % Liu, Y., X.S. Liang, and R.H. Weisberg, 2007: Rectification of the bias % in the wavelet power spectrum. Journal of Atmospheric and Oceanic % Technology, 24(12), 2093-2102. %
% Torrence, C., and G. P. Compo, 1998: A practical guide to wavelet % analysis. Bull. Amer. Meteor. Soc., 79, 61�78. %========================================================================
% Here starts with the original file header:

%WAVETEST Example Matlab script for WAVELET, using NINO3 SST dataset
%
% Written January 1998 by C. Torrence
%
% Modified Oct 1999, changed Global Wavelet Spectrum (GWS) to be sideways,
% changed all "log" to "log2", changed logarithmic axis on GWS
% a normal axis.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


Read Wavelet module This module is the translation from the matlab tools by Torrence and Compo (1998).

See the bottom of this notebook for the inside of this module.

In [2]:
from wavelet import wavelet


Construct sine waves
Section 5 in Liu et al. [2007]

In [3]:
dt=1./24. # 1 hour
t1=1.
t2=8.
t3=32.
t4=128.
t5=365.
n=365*24*13
t=np.arange(n)*dt
pi2=np.pi*2.

#5 sine waves
sst=np.sin(t/t1*pi2)+np.sin(t/t2*pi2)+np.sin(t/t3*pi2)+np.sin(t/t4*pi2)+np.sin(t/t5*pi2)

time = t/365. + 1993.0   # construct time array
xlim = [time[0],time[-1]]  # plotting range

plt.plot(time,sst)
xrange=plt.xlim(xlim)


Wavelet Spectrum by Torrence and Compo [1998]¶

Set wavelet parameters

In [4]:
pad = 1      # pad the time series with zeroes (recommended)
dj = 1./8.    # this will do 8 sub-octaves per octave
s0 = 6.*dt    # this says start at a scale of 6hour
j1 = 12./dj    # this says do 12 powers-of-two with dj sub-octaves each
mother = 'Morlet'


Wavelet transform

In [5]:
wave,period,scale,coi = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (np.abs(wave))**2         # compute wavelet power spectrum


Global wavelet spectrum

In [6]:
global_ws = power.sum(axis=1)/float(n)   # time-average over all times


Wavelet Spectrum by Liu et al [2007]¶

__Bias rectification__

divided by scales

In [7]:
########################
#  Spectrum
########################
powers=np.zeros_like(power)
for k in range(len(scale)):
powers[k,:] = power[k,:]/scale[k]

########################
#  Global Spectrum
########################
global_wss = global_ws/scale


PLOT¶

Corresponding to Fig. 2 in Liu et al. (2007)

In [10]:
#figure size
fig=plt.figure(figsize=(10,10))

# subplot positions
width= 0.65
height = 0.28;
pos1a = [0.1, 0.75, width, 0.2]
pos1b = [0.1, 0.4, width, height]
pos1c = [0.79, 0.4, 0.18, height]
pos1d = [0.1, 0.02, width, height]
pos1e = [0.79,0.02, 0.18, height]
#########################################
#---- a) Original signal
#########################################
ax.plot(time,sst,"r-")
ax.set_ylabel('')
plt.title('a) Original Sine curve')

#########################################
#   b) Wavelet spectrum
#########################################

#--- Contour plot wavelet power spectrum
Yticks = 2.**(np.arange(np.int(np.log2(np.min(period))),np.int(np.log2(np.max(period)))+1))
bx.contourf(time,np.log2(period),np.log2(power))
bx.set_xlabel('Time (year)')
bx.set_ylabel('Period (days)')
import matplotlib.ticker as ticker
ymajorLocator=ticker.FixedLocator(np.log2(Yticks))
bx.yaxis.set_major_locator( ymajorLocator )
ticks=bx.yaxis.set_ticklabels(Yticks)
plt.title('b) Wavelet Power Spectrum')

# cone-of-influence, anything "below" is dubious
ts = time;
coi_area = np.concatenate([[np.max(scale)], coi, [np.max(scale)],[np.max(scale)]])
ts_area = np.concatenate([[ts[0]], ts, [ts[-1]] ,[ts[0]]]);
L = bx.plot(ts_area,np.log2(coi_area),'k',linewidth=3)
F=bx.fill(ts_area,np.log2(coi_area),'k',alpha=0.3,hatch="x")

#########################################
#   c) Global Wavelet spectrum
#########################################

#--- Plot global wavelet spectrum
cx.plot(np.log2(global_ws),np.log2(period),"r-")
ylim=cx.set_ylim(np.log2([period.min(),period.max()]))
cx.invert_yaxis()
plt.title('c) Global Wavelet Spectrum')
xrangec=cx.set_xlim([0,1.25*np.max(np.log2(global_ws))])

#########################################
#   d) Wavelet spectrum
#########################################

#--- Contour plot wavelet power spectrum
dx.contourf(time,np.log2(period),np.log2(powers))
dx.set_xlabel('Time (year)')
dx.set_ylabel('Period (days)')
ymajorLocator=ticker.FixedLocator(np.log2(Yticks))
dx.yaxis.set_major_locator( ymajorLocator )
dx.yaxis.set_ticklabels(Yticks)
plt.title('d) Wavelet Power Spectrum')

# cone-of-influence, anything "below" is dubios
L = dx.plot(ts_area,np.log2(coi_area),'k',linewidth=3)
F=dx.fill(ts_area,np.log2(coi_area),'k',alpha=0.3,hatch="x")

#########################################
#   e) Global Wavelet spectrum
#########################################

#--- Plot global wavelet spectrum
ex.plot(np.log2(global_wss),np.log2(period),"r-")
ylim=ex.set_ylim(np.log2([period.min(),period.max()]))
ex.invert_yaxis()
plt.title('e) Global Wavelet Spectrum')
xrangec=ex.set_xlim([0,1.25*np.max(np.log2(global_wss))])

#########################################
#   Save fig
#########################################
plt.savefig("wavelet_test_sine.png")


List of Modules¶

In [11]:
!cat wavelet.py

def wavelet(Y,dt,pad=0.,dj=0.25,s0=-1,J1=-1,mother="MORLET",param=-1):
"""
This function is the translation of wavelet.m by Torrence and Compo

import wave_bases from wave_bases.py

The following is the original comment in wavelet.m

#WAVELET  1D Wavelet transform with optional singificance testing
%
%
%   Computes the wavelet transform of the vector Y (length N),
%   with sampling rate DT.
%
%   By default, the Morlet wavelet (k0=6) is used.
%   The wavelet basis is normalized to have total energy=1 at all scales.
%
%
% INPUTS:
%
%    Y = the time series of length N.
%    DT = amount of time between each Y value, i.e. the sampling time.
%
% OUTPUTS:
%
%    WAVE is the WAVELET transform of Y. This is a complex array
%    of dimensions (N,J1+1). FLOAT(WAVE) gives the WAVELET amplitude,
%    ATAN(IMAGINARY(WAVE),FLOAT(WAVE) gives the WAVELET phase.
%    The WAVELET power spectrum is ABS(WAVE)^2.
%    Its units are sigma^2 (the time series variance).
%
%
% OPTIONAL INPUTS:
%
% *** Note *** setting any of the following to -1 will cause the default
%               value to be used.
%
%    PAD = if set to 1 (default is 0), pad time series with enough zeroes to get
%         N up to the next higher power of 2. This prevents wraparound
%         from the end of the time series to the beginning, and also
%         speeds up the FFT's used to do the wavelet transform.
%         This will not eliminate all edge effects (see COI below).
%
%    DJ = the spacing between discrete scales. Default is 0.25.
%         A smaller # will give better scale resolution, but be slower to plot.
%
%    S0 = the smallest scale of the wavelet.  Default is 2*DT.
%
%    J1 = the # of scales minus one. Scales range from S0 up to S0*2^(J1*DJ),
%        to give a total of (J1+1) scales. Default is J1 = (LOG2(N DT/S0))/DJ.
%
%    MOTHER = the mother wavelet function.
%             The choices are 'MORLET', 'PAUL', or 'DOG'
%
%    PARAM = the mother wavelet parameter.
%            For 'MORLET' this is k0 (wavenumber), default is 6.
%            For 'PAUL' this is m (order), default is 4.
%            For 'DOG' this is m (m-th derivative), default is 2.
%
%
% OPTIONAL OUTPUTS:
%
%    PERIOD = the vector of "Fourier" periods (in time units) that corresponds
%           to the SCALEs.
%
%    SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J1
%            where J1+1 is the total # of scales.
%
%    COI = if specified, then return the Cone-of-Influence, which is a vector
%        of N points that contains the maximum period of useful information
%        at that particular time.
%        Periods greater than this are subject to edge effects.
%        This can be used to plot COI lines on a contour plot by doing:
%
%              contour(time,log(period),log(power))
%              plot(time,log(coi),'k')
%
%----------------------------------------------------------------------------
%   Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
%
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made. This
%   routine is provided as is without any express or implied warranties
%   whatsoever.
%
% Notice: Please acknowledge the use of the above software in any publications:
%    Wavelet software was provided by C. Torrence and G. Compo,
%      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
%
% Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
%            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
%
% Please send a copy of such publications to either C. Torrence or G. Compo:
%  Dr. Christopher Torrence               Dr. Gilbert P. Compo
%  Research Systems, Inc.                 Climate Diagnostics Center
%  4990 Pearl East Circle                 325 Broadway R/CDC1
%  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
%----------------------------------------------------------------------------"""
#modules
import numpy as np
from wave_bases import wave_bases

#set default
n1 = len(Y)
if (s0 == -1): s0=2.*dt
if (dj == -1): dj = 1./4.
if (J1 == -1): J1=np.fix((np.log(n1*dt/s0)/np.log(2))/dj)
if (mother == -1): mother = 'MORLET'
#print "s0=",s0
#print "J1=",J1

#....construct time series to analyze, pad if necessary
x = Y - np.mean(Y);
base2 = np.fix(np.log(n1)/np.log(2) + 0.4999)   # power of 2 nearest to N
temp=np.zeros(int((2**(base2+1)-n1),))
x=np.concatenate((x,temp))

n = len(x)

#....construct wavenumber array used in transform [Eqn(5)]
k = np.arange(1,np.fix(n/2)+1)
k = k*(2.*np.pi)/(n*dt)
k = np.concatenate((np.zeros((1,)),k, -k[-2::-1]));

#....compute FFT of the (padded) time series
f = np.fft.fft(x)    # [Eqn(3)]

#....construct SCALE array & empty PERIOD & WAVE arrays
scale=np.array([s0*2**(i*dj) for i in range(0,int(J1)+1)])
period = scale.copy()
wave = np.zeros((int(J1)+1,n),dtype=np.complex)  # define the wavelet array  # make it complex
# loop through all scales and compute transform
for a1 in range(0,int(J1)+1):
daughter,fourier_factor,coi,dofmin=wave_bases(mother,k,scale[a1],param)
wave[a1,:] = np.fft.ifft(f*daughter)  # wavelet transform[Eqn(4)]
period = fourier_factor*scale
coi=coi*dt*np.concatenate(([1.E-5],np.arange(1.,(n1+1.)/2.-1),np.flipud(np.arange(1,n1/2.)),[1.E-5])) # COI [Sec.3g]
wave = wave[:,:n1]  # get rid of padding before returning
return wave,period,scale,coi
# end of code

In [12]:
!cat wave_bases.py

def wave_bases(mother,k,scale,param):
"""
This is translation of wave_bases.m by Torrence and Gilbert P. Compo

The folloing is the original README

%    WAVE_BASES  1D Wavelet functions Morlet, Paul, or DOG
%
%  [DAUGHTER,FOURIER_FACTOR,COI,DOFMIN] = ...
%      wave_bases(MOTHER,K,SCALE,PARAM);
%
%   Computes the wavelet function as a function of Fourier frequency,
%   used for the wavelet transform in Fourier space.
%   (This program is called automatically by WAVELET)
%
% INPUTS:
%
%    MOTHER = a string, equal to 'MORLET' or 'PAUL' or 'DOG'
%    K = a vector, the Fourier frequencies at which to calculate the wavelet
%    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 of 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
%             (either 2 for Morlet and Paul, or 1 for the DOG)
%
%----------------------------------------------------------------------------
%   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
%   University of Colorado, Program in Atmospheric and Oceanic Sciences.
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.
%----------------------------------------------------------------------------
"""
#import modules
import numpy as np

#
mother = mother.upper()
n = len(k)
# define Heaviside step function
def ksign(x):
y=np.zeros_like(x)
y[x>0]=1
return y
#
if mother=='MORLET':  #-----------------------------------  Morlet
if (param == -1): param = 6.
k0 = param
expnt = -(scale*k - k0)**2/2. *ksign(k)
norm = np.sqrt(scale*k[1])*(np.pi**(-0.25))*np.sqrt(n)    # total energy=N   [Eqn(7)]
daughter = norm*np.exp(expnt)
daughter = daughter*ksign(k)  # Heaviside step function
fourier_factor = (4.*np.pi)/(k0 + np.sqrt(2. + k0**2)) # Scale-->Fourier [Sec.3h]
coi = fourier_factor/np.sqrt(2)            # Cone-of-influence [Sec.3g]
dofmin = 2.                          # Degrees of freedom
elif mother=='PAUL': #--------------------------------  Paul
if (param == -1): param = 4.
m = param
expnt = -(scale*k)*ksign(k)
norm = np.sqrt(scale*k[1])*(2.**m/np.sqrt(m*np.prod(np.arange(2,2*m))))*np.sqrt(n)
daughter = norm*((scale*k)**m)*np.exp(expnt)
daughter = daughter*ksign(k)      # Heaviside step function
fourier_factor = 4*np.pi/(2.*m+1.)
coi = fourier_factor*np.sqrt(2)
dofmin = 2.
elif mother=='DOG':  #--------------------------------  DOG
if (param == -1): param = 2.
m = param
expnt = -(scale*k)**2 / 2.0
from scipy.special import gamma
norm = np.sqrt(scale*k[1]/gamma(m+0.5))*np.sqrt(n)
daughter = -norm*(1j**m)*((scale*k)**m)*np.exp(expnt);
fourier_factor = 2.*np.pi*np.sqrt(2./(2.*m+1.))
coi = fourier_factor/np.sqrt(2)
dofmin = 1.
else:
raise Exception("Mother must be one of MORLET,PAUL,DOG")

return daughter,fourier_factor,coi,dofmin

# end of code

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: