%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.load 'sst_nino3.dat'   % input SST time seriessst = 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 arrayxlim = [1870,2000];  % plotting rangepad = 1;      % pad the time series with zeroes (recommended)dj = 0.25;    % this will do 4 sub-octaves per octaves0 = 2*dt;    % this says start at a scale of 6 monthsj1 = 7/dj;    % this says do 7 powers-of-two with dj sub-octaves eachlag1 = 0.72;  % lag-1 autocorrelation for red noise backgroundmother = '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) arraysig95 = power ./ sig95;         % where ratio > 1, power is significant% Global wavelet spectrum & significance levels:global_ws = variance*(sum(power')/n);   % time-average over all timesdof = n - scale;  % the -scale corrects for padding at edgesglobal_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother);% Scale-average between El Nino periods of 2--8 yearsavg = find((scale >= 2) & (scale < 8));Cdelta = 0.776;   % this is for the MORLET waveletscale_avg = (scale')*(ones(1,n));  % expand scale --> (J+1)x(N) arrayscale_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%------------------------------------------------------ Plotting%--- Plot time seriessubplot('position',[0.1 0.75 0.65 0.2])plot(time,sst)set(gca,'XLim',xlim(:))xlabel('Time (year)')ylabel('NINO3 SST (degC)')title('a) NINO3 Sea Surface Temperature (seasonal)')hold off%--- Contour plot wavelet power spectrumsubplot('position',[0.1 0.37 0.65 0.28])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));  %*** or use 'contourfill'%imagesc(time,log2(period),log2(power));  %*** uncomment for 'image' plotxlabel('Time (year)')ylabel('Period (years)')title('b) NINO3 SST Wavelet Power Spectrum')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 oncontour(time,log2(period),sig95,[-99,1],'k');hold on% cone-of-influence, anything "below" is dubiousplot(time,log2(coi),'k')hold off%--- Plot global wavelet spectrumsubplot('position',[0.77 0.37 0.2 0.28])plot(global_ws,log2(period))hold onplot(global_signif,log2(period),'--')hold offxlabel('Power (degC^2)')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_ws)])%--- Plot 2--8 yr scale-average time seriessubplot('position',[0.1 0.07 0.65 0.2])plot(time,scale_avg)set(gca,'XLim',xlim(:))xlabel('Time (year)')ylabel('Avg variance (degC^2)')title('d) 2-8 yr Scale-average Time Series')hold onplot(xlim,scaleavg_signif+[0,0],'--')hold off% end of code