% 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 to % a normal axis. clear all; close all; 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 % Significance levels: (variance=1 for the normalized SST) [signif,fft_theor] = wave_signif(1.0,dt,scale,0,lag1,-1,-1,mother); sig95 = (signif')*(ones(1,n)); % expand signif --> (J+1)x(N) array sig95 = power ./ sig95; % where ratio > 1, power is significant % Global wavelet spectrum & significance levels: global_ws = variance*(sum(power')/n); % time-average over all times dof = n - scale; % the -scale corrects for padding at edges global_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother); % Scale-average between El Nino periods of 2--8 years avg = find((scale >= 2) & (scale < 8)); Cdelta = 0.776; % this is for the MORLET wavelet scale_avg = (scale')*(ones(1,n)); % expand scale --> (J+1)x(N) array scale_avg = power ./ scale_avg; % [Eqn(24)] scale_avg = variance*dj*dt/Cdelta*sum(scale_avg(avg,:)); % [Eqn(24)] scaleavg_signif = wave_signif(variance,dt,scale,2,lag1,-1,[2,7.9],mother); % whos %------------------------------------------------------ Plot: hight = 0.2; pos1a = [0.1 0.37 0.65 hight]; pos1b = [0.78 0.37 0.15 hight]; pos2a = [0.1 0.14 0.65 hight]; pos2b = [0.78 0.14 0.15 hight]; figure('visible','off'); orient tall; %--- Contour plot wavelet power spectrum subplot('position',pos1a) levels = [0.0625,0.125,0.25,0.5,1,2,4,8,16] ; Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period)))); contour(time,log2(period),log2(power),log2(levels)); caxis([log2(levels(1)), log2(levels(end))]); % xlabel('Time (year)') ylabel('Period (years)') set(gca,'XLim',xlim(:)) set(gca,'YLim',log2([min(period),max(period)]), ... 'YDir','reverse', ... 'YTick',log2(Yticks(:)), ... 'YTickLabel',Yticks) % 95% significance contour, levels at -99 (fake) and 1 (95% signif) hold on [cs,h] = contour(time,log2(period),sig95,[-99,1],'k'); set(h,'linewidth',1); hold on; % cone-of-influence, anything "below" is dubious % plot(time,log2(coi),'k'); hold off ts = time; coi_area = [max(scale) coi max(scale) max(scale)]; ts_area = [ts(1) ts ts(length(ts)) ts(1)]; L = plot(ts_area,log2(coi_area),'k'); set(L,'linewidth',.3); hold on hatch(L,45,'k','-',5,.3); hold on hatch(L,135,'k','-',5,.3); hold on %--- Plot global wavelet spectrum subplot('position',pos1b) plot(global_ws,log2(period)); hold on; % plot(global_signif,log2(period),'--'); hold off % xlabel('Power (degC^2)') set(gca,'YLim',log2([min(period),max(period)]), ... 'YDir','reverse', ... 'YTick',log2(Yticks(:)), ... 'YTickLabel','') set(gca,'XLim',[0,1.25*max(global_ws)]) %================== Bias rectification start ========================= %--- divided by scales: for k=1:length(scale) powers(k,:) = power(k,:)/scale(k); end global_wss = global_ws./scale; c = powers; c = c(:); %================== Bias rectification end ========================= %--- Contour plot wavelet power spectrum subplot('position',pos2a) levelss = 2.^[-6:2] ; Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period)))); contour(time,log2(period),log2(powers),log2(levelss)); caxis([log2(levelss(1)), log2(levelss(end))]); xlabel('Time (year)') ylabel('Period (years)') set(gca,'XLim',xlim(:)) set(gca,'YLim',log2([min(period),max(period)]), ... 'YDir','reverse', ... 'YTick',log2(Yticks(:)), ... 'YTickLabel',Yticks) % 95% significance contour, levels at -99 (fake) and 1 (95% signif) hold on [cs,h] = contour(time,log2(period),sig95,[-99,1],'k'); hold on set(h,'linewidth',1); hold on; % cone-of-influence, anything "below" is dubious % plot(time,log2(coi),'k'); hold off L = plot(ts_area,log2(coi_area),'k'); set(L,'linewidth',.3); hold on hatch(L,45,'k','-',5,.3); hold on hatch(L,135,'k','-',5,.3); hold on %--- Plot global wavelet spectrum subplot('position',pos2b) plot(global_wss,log2(period)) hold on % title('c) Global Wavelet Spectrum') set(gca,'YLim',log2([min(period),max(period)]), ... 'YDir','reverse', ... 'YTick',log2(Yticks(:)), ... 'YTickLabel','') set(gca,'XLim',[0,1.25*max(global_wss)]) %% xlabel('Global') xlabel('Time Averaged') print -djpeg wavelet_test_ElNino3.jpg; close all;