- 积分
 - 11
 
	- 贡献
 -  
 
	- 精华
 
	- 在线时间
 -  小时
 
	- 注册时间
 - 2014-3-10
 
	- 最后登录
 - 1970-1-1
 
 
 
 
 
 
 | 
	
 
 
 楼主 |
发表于 2014-5-5 21:25:47
|
显示全部楼层
 
 
 
 
st = [ 
 (此处为数据) 
] 
variance = std(sst)^2; 
sst = (sst - mean(sst))/sqrt(variance) ; 
 
n = length(sst); 
dt = 1 ; 
time = [0:length(sst)-1]*1 + 1951.0 ;  % construct time array 
xlim = [1950,2005];  % plotting range 
pad = 1;      % pad the time series with zeroes (recommended) 
dj = 0.25;    % this will do 4 sub-octaves per octave 
s0 = 1*dt;    % this says start at a scale of 6 months 
j1 = 6/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); 
 
%------------------------------------------------------ Plotting 
 
%--- Plot time series 
subplot('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 spectrum 
subplot('position',[0.1 0.30 0.45 0.55]) 
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' plot 
xlabel('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 on 
contour(time,log2(period),sig95,[-99,1],'k'); 
hold on 
% cone-of-influence, anything "below" is dubious 
plot(time,log2(coi),'k') 
hold off 
 
%--- Plot global wavelet spectrum 
subplot('position',[0.57 0.30 0.1 0.55]) 
plot(global_ws,log2(period)) 
hold on 
plot(global_signif,log2(period),'--') 
hold off 
xlabel('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 series 
subplot('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 on 
plot(xlim,scaleavg_signif+[0,0],'--') |   
 
 
 
 |