sst = nianpj; % precipitation为保存在MATLAB workspace中的年降水量矢量数据
variance = std(sst)^2;%计算方差
sst = (sst - mean(sst))/sqrt(variance) %数据标准化。如果导入的降水量数据已经标准化，此步骤可忽略
n = length(sst); %降水量数据的矢量长度
dt = 1%设置时间步长为1年
time = [0:n-1]*dt+ 1981 ;%设置时间矢量，从1961到
xlim = [1981,];%绘图时x轴范围
pad = 1; %填补模式为1。本时间序列长度为47，不足26，故选择0填补
dj = 0.125; % this will do 4 sub-octaves per octave
s0 = 2*dt; %可分辨的最小尺度为2年
j1 = 6/dj; % 为26中的6，即以2的幂指数，根据时间序列长度确定
lag1 = 0.72; %红噪声背景谱的自相关检验参数
mother = 'Morlet';%母小波为“Morlet”
%%-------------小波转换-------------------------------------------
%计算小波系数、周期、尺度，请参考软件中该函数说明
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (abs(wave)).^2% 计算小波功率谱
%%-------显著性检验 (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; %当比值大于1，表明强度值显著
%%------------------小波全谱和显著性检验--------------------------
global_ws = variance*(sum(power')/n); %所有时间尺度上的平均
dof = n - scale; % 边缘效应的尺度订正
global_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother);
%%-------------保存小波分析结果--------------------------
save period.dat period -ascii;
save scale.dat scale -ascii;
save power.dat power -ascii;
save global_ws.dat global_ws -ascii;
save global_signif.dat global_signif -ascii;
%%------绘图----------------------------------------------
%--- 绘制时间序列图------------
subplot('position',[0.1 0.75 0.65 0.2])
plot(time,sst)
set(gca,'XLim',xlim(:))
xlabel('Time (day)')
ylabel('S(l/s)')
title('a) 标准化后的年降水时序 (seasonal)')
hold off
%--- 绘制小波功率谱等值线图----------------------
subplot('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' plot
xlabel('年份')
ylabel('周期 (年)')
title('b) 年降水小波功率谱')
set(gca,'XLim',xlim(:))
set(gca,'YLim',log2([min(period),max(period)]), ...
'YDir','reverse', ...
'YTick',log2(Yticks(:)), ...
'YTickLabel',Yticks)
% Contourfill
imagesc(time,log2(period),log2(power));
xlabel('年份')
ylabel('周期 (年)')
title('b) 年降水小波功率谱')
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
% 影响锥线，锥线一下的周期特征存在不确定性
plot(time,log2(coi),'k')
hold off
%--- 绘制小波全谱---------------
subplot('position',[0.77 0.37 0.2 0.28])
plot(global_ws,log2(period))
hold on
plot(global_signif,log2(period),'--')
hold off
xlabel('功率 (l/s^2)')
title('c) 小波全谱')
set(gca,'YLim',log2([min(period),max(period)]), ...
'YDir','reverse', ...
'YTick',log2(Yticks(:)), ...
'YTickLabel','')
set(gca,'XLim',[0,1.25*max(global_ws)])
%---补充绘制小波全谱图(反向坐标轴)----------------------
subplot('position',[0.1 0.1 0.2 0.15])
plot(log2(period),global_ws)
hold on
plot(log2(period),global_signif,'--')
hold off
xlabel('周期 (年)')
ylabel('功率(1/s^2)')
title('d) 小波全谱')
set(gca,'XLim',log2([min(period),max(period)]),'XTick',log2(Yticks(:)), ...
'XTickLabel',Yticks)
以上源代码运行后提示
Undefined function 'wavelet' for input arguments of type 'double'.
Error in xioabofenxi (line 19)
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
求助高手给我提供一个解决方法
