- 积分
- 709
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-11-9
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
最近尝试做小波分析,在论坛里找到几个程序,尝试照葫芦画瓢,但是总有很多不明白的地方。
比如:小波系数实部等值线图和功率谱图很相似?方差和总功率又是什么关系?
总之现在是很困惑,夜里睡不着出来做个东西,总是不能和例子一样,很郁闷。
目前最想请教高人指点一下下面的图是怎么回事:
这是我根据几个程序自己写的,原代码如下:
clear
load('xxx.txt')
time=xxx(:,1);
data=xxx(:,2);
xx=gaussfilter(data,9);
sst =xx;
n=length(xx);
variance = std(sst)^2;
sst = (sst - mean(sst))/sqrt(variance);
dt = 1;
year = [0:n-1]+ 1900.0 ;
xlim = [1900,2009];
pad = 1;
dj = 1/12;
s0 = 10*dt;
j1 = 6.5/dj;
lag1 = 0.72;
mother = 'Morlet';
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (abs(wave)).^2 ;
global_ws =variance*(sum(power')/n);
[signif,fft_theor] = wave_signif(1.0,dt,scale,0,lag1,-1,-1,mother);
sig95 = (signif')*(ones(1,n));
sig95 = power ./ sig95;
dof = n - scale;
global_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother);
figure(1)
subplot(2,1,1)
levels = [0,0.5,1.0,1.5,2.0,2.5];
v = [0,0.5,1.0,1.5];
Yticks = 10:10:100;
[c1,h1]=contourf(year,period,real(wave),levels,'k-');
clabel(c1,h1,v,'fontsize',5);
xlabel('年份/year')
ylabel('周期/年 period/year')
title('小波功率谱')
set(gca,'XLim',xlim(:))
set(gca,'YLim',[10 100], ...
'YDir','default', ...
'YTick',Yticks(:), ...
'YTickLabel',Yticks)
hold on
levels = [-0.5,-1.0,-1.5,-2.0,-2.5];
v = [-0.5,-1.0,-1.5];
[c,h] = contour(year,period,power,levels,'r--');
clabel(c,h,v,'fontsize',5);
contour(year,period,sig95,[-99,1],'g');
plot(year,coi,'g')
hold off
subplot(2,1,2)
plot(period,global_ws)
levels= [10,20,30,40,50,60,70,80,90,100];
hold on
plot(period,global_signif,'--')
title('Global Wavelet Spectrum')
set(gca,'XLim',[10,100], ...
'XTick',levels,...
'XTickLabel',levels)
xlabel('Time (year)')
ylabel('Power ')
hold off
也许这只是很低级的问题,让大家见笑了,一时想不同,找不到原因,还望高人赐教。
|
|