Continuous Wavelet Spectrum

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
%
% See "http://paos.colorado.edu/research/wavelets/"
% 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.

Reading modules and data

Read usual modules

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).
http://paos.colorado.edu/research/wavelets/software.html

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 [8]:
#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=fig.add_axes(pos1a)
ax.plot(time,sst,"r-")
ax.set_ylabel('')
plt.title('a) Original Sine curve')

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

#--- Contour plot wavelet power spectrum
bx=fig.add_axes(pos1b,sharex=ax) 
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=fig.add_axes(pos1c,sharey=bx)
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=fig.add_axes(pos1d,sharex=ax)
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=fig.add_axes(pos1e,sharey=dx)
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")
-c:31: RuntimeWarning: divide by zero encountered in log2
-c:65: RuntimeWarning: divide by zero encountered in log2

List of Modules

In [9]:
!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
%
%   [WAVE,PERIOD,SCALE,COI] = wavelet(Y,DT,PAD,DJ,S0,J1,MOTHER,PARAM)
%
%   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
%  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
%----------------------------------------------------------------------------"""  
    #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);
    if (pad == 1):
        base2 = np.fix(np.log(n1)/np.log(2) + 0.4999)   # power of 2 nearest to N
        temp=np.zeros((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 [10]:
!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