clear all
load 'sst.txt'
% 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 =1 ;
time = [0:length(sst)-1]dt + 1971.0 ; % construct time array
xlim = [1971,2020]; % plotting range
pad = 1; % pad the time series with zeroes (recommended)
dj = 0.25; % this will do 4 sub-octaves per octave
s0 = 2dt; % this says start at a scale of 24 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

modulus=abs(wave); %計(jì)算小波系數(shù)的模
variance1=sum(power')/n;%計(jì)算小波方差
%畫小波系數(shù)實(shí)部等值線圖
% subplot(3,1,1)
levels = [-4,-3.5,-3,-2.5,-2.0,-1.5,-1.0,-0.5,0,0.5,1.0,1.5,2.0,2.5,3,3.5,4];
v = [-3.5,-3,-2.5,-2.0,-1.5,-1.0,-0.5,0,0.5,1.0,1.5,2.0,2.5,3,3.5];
Yticks = 0:10:50;
[c,h]=contourf(time,period,real(wave),levels,'k-');
clabel(c,h,v,'fontsize',5);
y=real(wave)
xlabel('年份/year')
ylabel('周期/年 period/year')
set(gca,'XLim',xlim(:))
set(gca,'YLim',[0 50], ...
'YDir','default', ...
'YTick',Yticks(:), ...
'YTickLabel',Yticks)
hold on
% 畫小波方差圖
subplot(3,1,2)
plot(period,variance1,'k-')
hold on;
levels= [0,10,20,30,40,50];
title('(b)')
set(gca,'XLim',[1 50], ...
'XTick',levels,...
'XTickLabel',levels)
xlabel('周期/a')
ylabel('方差 variance')
hold on

% 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 = variancedjdt/Cdeltasum(scale_avg(avg,:)); % [Eqn(24)]
scaleavg_signif = wave_signif(variance,dt,scale,2,lag1,-1,[2,7.9],mother);

whos
% period=period(1:end-2);
% Yticks=Yticks(1:end-2);
%------------------------------------------------------ Plotting
levels = [0.0625,0.125,0.25,0.5,1,2,4,8,16] ;
Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
load('1.mat');
colormap(map(140:240,:))
contourf(time,log2(period),log2(power),log2(levels)); %*** or use 'contourfill'
%imagesc(time,log2(period),log2(power)); %*** uncomment for 'image' plot
xlabel('Time (year)','FontSize',14)
ylabel('Period (years)','FontSize',14)
title(' Wavelet Power Spectrum','FontSize',14)
set(gca,'XLim',[1959 2018]);
period=sort(period);
set(gca,'YLim',log2([min(period),max(period)]), ...
'YDir','reverse', ...
'YTick',log2(Yticks(:)), ...
'YTickLabel',Yticks)
hc = colorbar;
set(hc,'fontweight','Bold')
set(gca,'fontweight','Bold','FontSize',14);
% 95% significance contour, levels at -99 (fake) and 1 (95% signif)
hold on
contour(time,log2(period),sig95,[-99,1],'b','linewidth',2);
hold on
% cone-of-influence, anything "below" is dubious
plot(time,log2(coi),'k','linewidth',2)
hold off
axis xy
set(gca,'linewi',1)
set(gca,'YLim',[1 5])