First written on Nov,19, 2014:
Author: tmiyama at gmail
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:
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. http://ocgweb.marine.usf.edu/~liu/wavelet.html
Torrence, C., and G. P. Compo, 1998: A practical guide to wavelet analysis. Bull. Amer. Meteor. Soc., 79, 6178. http://paos.colorado.edu/research/wavelets/
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: yliu18@gmail.com
% 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.
Read usual modules
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.
from wavelet import wavelet
Construct sine waves
Section 5 in Liu et al. [2007]
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)
Set wavelet parameters
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
wave,period,scale,coi = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (np.abs(wave))**2 # compute wavelet power spectrum
Global wavelet spectrum
global_ws = power.sum(axis=1)/float(n) # time-average over all times
Bias rectification
divided by scales
########################
# 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
Corresponding to Fig. 2 in Liu et al. (2007)
#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")
!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(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
!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