爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3858|回复: 4

向大神请教小波分析的问题

[复制链接]

新浪微博达人勋

发表于 2019-7-18 16:11:09 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
向大神请教:做了5年,每日数据的小波分析,但是Scale(纵坐标)最大是256,怎么设置成512?

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

% 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);


%--- Contour plot wavelet power spectrum
%Yticks1=[4,8,16,32,64,128,256,512]
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
XTick1 = [2013,2014,2015,2016,2017]
xlabel('Time (year)')
ylabel('Period (days)')
%colorbar('EASTOUTSIDE')
title('ACI Wavelet Power Spectrum')
%set(gca,'xticklabel',{'2015/01','2015/07', '2016/01', '2016/07','2017/01', '2017/07'})      
%set(gca,'XLim',xlim(:),...
%    'XTick',xlim/365+2015, ...
%        'XTickLabel',XTick1)
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
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-7-18 16:14:26 | 显示全部楼层
本帖最后由 ellen 于 2019-7-18 16:15 编辑

图片是这样的
1563334213(1).png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-19 08:20:10 | 显示全部楼层
改一下Ylim。话说lag1自相关那里也要改,不然显著性检验有问题。你去好好看看Torrence1998BAMS那篇文章吧
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-19 10:24:26 | 显示全部楼层
请问您数据输入的格式是什么样子的,我输入的始终都有问题,可以传一份给我吗…… 谢谢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-19 11:20:14 | 显示全部楼层
j1 = 7/dj; 7改成8就行了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表