Replicate the example of Torrence & Compo 1998

[pinshuai]

Hi, I am trying to replicate the results from the classic wavelet example of Torrence & Compo, 1998 using pywt.cwt(), but could not get it right. I have pasted the original figure produced using Matlab in Torrence & Compo's paper. Their source codes can be found here. A code snippet is copied here:

load 'sst_nino3.dat'   % input SST time series
sst = sst_nino3;

%------------------------------------------------------ Computation

% normalize by standard deviation (not necessary, but makes it easier
% to compare with plot on Interactive Wavelet page, at
% "http://paos.colorado.edu/research/wavelets/plot/"
variance = std(sst)^2;
sst = (sst - mean(sst))/sqrt(variance) ;

n = length(sst);
dt = 0.25 ;
time = [0:length(sst)-1]*dt + 1871.0 ;  % construct time array
xlim = [1870,2000];  % plotting range
pad = 1;      % pad the time series with zeroes (recommended)
dj = 0.25;    % this will do 4 sub-octaves per octave
s0 = 2*dt;    % this says start at a scale of 6 months
j1 = 7/dj;    % this says do 7 powers-of-two with dj sub-octaves each
lag1 = 0.72;  % lag-1 autocorrelation for red noise background
mother = 'Morlet';

% Wavelet transform:
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (abs(wave)).^2 ;        % compute wavelet power spectrum

I am aware of an example code using pywt for the same dataset (https://github.com/PyWavelets/pywt/blob/master/demo/cwt_analysis.py). But it was using a different wavelet function (i.e., cmor) instead of the original Morlet function. Also the scales used are not the same. Below is the figure from the example code in pywt.

url = 'http://paos.colorado.edu/research/wavelets/wave_idl/nino3sst.txt'
sst = np.genfromtxt(url, skip_header=19)

t0 = 1871.0
dt = 0.25  # In years

# We also create a time array in years.
N = sst.size
time = np.arange(0, N) * dt + t0

# time, sst = pywt.data.nino()
# dt = time[1] - time[0]

# Taken from http://nicolasfauchereau.github.io/climatecode/posts/wavelet-analysis-in-python/
wavelet = 'cmor1.5-1.0'

scales = np.arange(1, 128)

[cfs, frequencies] = pywt.cwt(sst, scales, wavelet, dt)
power = (abs(cfs)) ** 2

period = 1. / frequencies
levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8]
f, ax = plt.subplots(figsize=(15, 10))
ax.contourf(time, np.log2(period), np.log2(power), np.log2(levels),
            extend='both')

ax.set_title('%s Wavelet Power Spectrum (%s)' % ('Nino1+2', wavelet))
ax.set_ylabel('Period (years)')
Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),
                        np.ceil(np.log2(period.max())))
ax.set_yticks(np.log2(Yticks))
ax.set_yticklabels(Yticks)
ax.invert_yaxis()
ylim = ax.get_ylim()
ax.set_ylim(ylim[0], -1)

If I change the wavelet function to morlet (wavelet = 'morl'), I got the following plot, which looks quite different from the original plot.

So, is cmor1.5-1.0 actually similar to Morlet used in the Matlab code? How to choose the scales in pywt so that I can replicate the Matlab example?

[grlee77]

Taking a quick look at the following section of Torrence's code: https://github.com/chris-torrence/wavelets/blob/master/wave_python/waveletFunctions.py#L170-L182

It appears that he defines the wavelet directly in the frequency domain (whereas the corresponding formulas in PyWavelets are in the time domain). The fact that the waveform is being multiplied by a Heaviside (step) function in frequency is going to make it complex-valued in the time domain (If you recall FFT theory, real-valued functions would have conjugate-symmetry in the frequency domain)

Based on that, I would say that "Morlet" in Torrence's toolbox corresponds to "cmor" in PyWavelets. I haven't tried to verify if the bandwidth / center frequency constants exactly match cmor1.5-1.0, but it is definitely a "cmor" wavelet and not a "morl" wavelet.

[pinshuai]

Thank you for the clarification. I am not very familiar with the FFT theory, but I think they must use a cmor wavelet since the plot generated by cmor1.5-1.0 look very similar to the original spectral plot.

In the original Matlab code, they used this equation to generate the scales.

In your example, you used scales = np.arange(1, 128) to generate a uniform distributed scales. Do you have any reasoning behind that? Or do you have any suggestions on choosing appropriate scales? Thanks.

[khider]

Having the same issue. I'm not sure how to choose: 1. the scales based on a frequency vector, 2. how to align the two outputs. Any help on this would be appreciated.