爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
123
返回列表 发新帖
楼主: 看云淡风清

[讨论] matlab小波分析图比较

[复制链接]

新浪微博达人勋

 楼主| 发表于 2017-12-25 15:28:14 | 显示全部楼层
晴晴 发表于 2017-12-25 12:29
还有我想请问一下我用自己数据为什么没有图一上面的图,下面的图2也没有颜色,好像用自带的数据就可以画出 ...

直接贴程序吧,你这么说我也不知道什么问题
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-26 09:30:09 | 显示全部楼层
看云淡风清 发表于 2017-12-25 15:28
直接贴程序吧,你这么说我也不知道什么问题

%WAVETEST Example Matlab script for WAVELET, using NINO3 SST dataset
%
% See "http://paos.colorado.edu/research/wavelets/"
% Written January 1998 by C. Torrence
%
% Modified Oct 1999, changed Global Wavelet Spectrum (GWS) to be sideways,
%   changed all "log" to "log2", changed logarithmic axis on GWS to
%   a normal axis.

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 = 0.25 ;
time = [0:length(sst)-1]*dt + 1871.0 ;  % construct time array
xlim = [1870,2000];  % 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);

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

whos

%------------------------------------------------------ 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.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('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.77 0.37 0.2 0.28])
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],'--')
hold off

% end of code
你好这个就是那个wavetest程序,就是我从家园下载下来的,我画出来的图和您的不一样,一个是您最上面的两个图没有,第三个图也没有颜色,不知道我是不是下错程序了,或者您方便给我一个您的程序,或者咱们加一个qq或者微信,1004368879@qq。com,也是微信号
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-12-26 17:05:50 | 显示全部楼层
晴晴 发表于 2017-12-26 09:30
%WAVETEST Example Matlab script for WAVELET, using NINO3 SST dataset
%
% See "http://paos.colora ...

两个图不是一个程序画的,一个是小波实部等值线和方差图,一个是功率谱图,上面这个程序画的是功率谱图,论坛有很多小波的程序,多看看就明白了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-27 10:45:29 | 显示全部楼层
看云淡风清 发表于 2017-12-26 17:05
两个图不是一个程序画的,一个是小波实部等值线和方差图,一个是功率谱图,上面这个程序画的是功率谱图, ...

请问一下,我看您的功率谱图有颜色是在哪个地方改的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-12-27 19:05:27 | 显示全部楼层
晴晴 发表于 2017-12-27 10:45
请问一下,我看您的功率谱图有颜色是在哪个地方改的

matlab可以设置颜色啊
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-27 22:49:47 | 显示全部楼层
学到了学到了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-3-7 19:59:56 | 显示全部楼层
谢谢,最近也需要画这样的图
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-3-28 23:42:04 | 显示全部楼层
请教楼主:功率谱图的纵坐标怎么设置的?我的数据系列是1966-2004年
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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