比较全的小波分析，小波系数、小波方差、模、模方，主周期，能量谱

vihenry2004

半年前开始学习小波分析，资料都是在网上找的，主要是气象家园的帖子，我把收集到的小波系数、小波方差、模、模方，主周期，能量谱（信度检验）都画了出来。分享出来跟大家交流。主要目的是感谢气象家园提供的学习平台，本版关于小波分析的帖子我基本都看过了，在此不一一列出，向所有帖主和参与讨论的网友致谢！我做的只是把我能理解的东西综合到了一起，半年前我还对小波分析一无所知，现在也有很多地方没有弄明白，所以以下内容里肯定有很多问题或错误，希望路过的高人不吝赐教！这也是本帖的目的之一。

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

以下是正文

数据来源：https://wenku.baidu.com/view/226 ... c4bb4cf7ecb722.html【这个例子的好处是有数据，有分析过程，还有最终画出的图，很全面】

一、数据处理，减小边界效应【网上看到的大部分程序没有这个步骤，上来就是小波变换，导致画出来的图与例子中的差别很大】

1.将原始数据按反年代-年代-反年代排列，例如：2005，2004，……，1966，1966，……2004，2005，2005，2004，……，1966。将排列好的数据复制到“q.txt”文件中，放在桌面。

2.打开matlab程序，在“命令窗口”里输入“loadC:\Users\lenovo\Desktop\q.txt”，回车，左上栏workspace出现黄色文件，将q另保存为“q.mat“。

3.在“命令窗口”里输入wavemenu，选择“One-Dimensional”下的子菜单“Complex Continuous Wavelet 1-D”，打开一维复连续小波界面，单击“File”菜单下的“Load Signal”按钮，载入径流时间序列q.mat。图的左侧为信号显示区域，右侧区域给出了信号序列和复小波变换的有关信息和参数，主要包括数据长度（Data Size）、小波函数类型（Wavelet：cgau、shan、fbsp和cmor）、取样周期（Sampling Period）、周期设置（Scale Setting）和运行按钮（Analyze），以及显示区域的相关显示设置按钮。本例中，我们选择cmor (1-1.5)、取样周期（Sampling Period）为1、最大尺度为32，单击“Analyze”运行按钮，计算小波系数。

4.然后单击“File”菜单下的“Save Coefficients”，保存为“confs.mat”。

5. matlab-file，打开“confs.mat“，左上角workspace出现四个文件，双击confs打开（总共排列年代不能大于1024年，即1024/3=341年，超过341年则要对原始排列序列进行截取），将前后39列删除，剩下中间原始的39列【根据自己的数据个数筛选，本例子是39年的数据】。

6. 另存为eg.mat

7. 为了绘制能谱，进行信度检验，另将没有处理过的原数据存为qq.txt（只要径流数据列，不要年份）

二、代码部分

clc;

clear;

dt=1;

t=[0:39-1]*dt+1966.0;%根据自己的数据修改

%进行连续小波变换得到小波系数矩阵,选择复morlet小波函数

load eg.mat;%载入小波系数

wf=coefs;

% 求得系数的实部

shibu=real(wf);

subplot(221);

contourf(shibu,20,'-');%中间的数值可调，值越大细节越详细

colorbar;

time=1970:5:2000;  % 按照自己的时间序列进行修改

xlabel('年份');

ylabel('时间尺度/年');

set(gca,'XTickLabel',time) %更新XTickLabel

% 小波方差是模的平方的算数平均

mo=abs(wf);% 求模

mofang=mo.^2;% 模的平方

fangcha=sum(mofang')/39;%小波方差

subplot(222);

plot(fangcha,'k-','linewidth',1.5);

xlabel('时间尺度/年');

ylabel('小波方差');

subplot(223);

contourf(mo,10,'-');%中间的数值可调，值越大细节越详细

time=1970:5:2000;  % 按照自己的时间序列进行修改

xlabel('年份');

ylabel('模');

colorbar;

set(gca,'XTickLabel',time) %更新XTickLabel

subplot(224);

contourf(mofang,10,'-');%中间的数值可调，值越大细节越详细

time=1970:5:2000;  % 按照自己的时间序列进行修改

xlabel('年份');

ylabel('模方');

colorbar;

set(gca,'XTickLabel',time) %更新XTickLabel

figure(2);% 小波实部过程线

subplot(221);

plot(t,shibu(28,:),'k-','linewidth',1.5);%第一主周期

hold on

plot([1965,2005],[0,0],'k-','linewidth',0.2);%y=0直线，按照自己的时间序列进行修改

ylabel('小波系数');

legend('28a特征时间尺度');

subplot(222);

plot(t,shibu(14,:),'k-','linewidth',1.5);%第二主周期

hold on

plot([1965,2005],[0,0],'k-','linewidth',0.2);%y=0直线;

ylabel('小波系数');

legend('14a特征时间尺度');

subplot(223);

plot(t,shibu(4,:),'k-','linewidth',1.5);%第三主周期

hold on

plot([1965,2005],[0,0],'k-','linewidth',0.2);%y=0直线

xlabel('年份');

ylabel('小波系数');

legend('4a特征时间尺度');

% 绘制能谱检验【这部分不懂，考别人的程序，不会改；没有用减小边界效应的eg.mat；期待高人指点】

load 'qq.txt'   % 载入数据序列

sst = qq;

variance =std(sst)^2;

sst = (sst -mean(sst))/sqrt(variance) ;

n = length(sst);

dt = 1 ;

time = [0:n-1]*dt +1966;

xlim = [1966,2004];

pad = 1;     

dj = 0.125;   

s0 = 2*dt;  

j1 = 6/dj;   

lag1 = 0.72;

mother = 'Morlet';

[wave,period,scale,coi]= wavelet(sst,dt,pad,dj,s0,j1,mother);

power =(abs(wave)).^2 ;  

[signif,fft_theor]= wave_signif(1.0,dt,scale,0,lag1,-1,-1,mother);

sig95 =(signif')*(ones(1,n));

sig95 = power ./sig95;   

global_ws =variance*(sum(power')/n);

dof = n - scale;

global_signif =wave_signif(variance,dt,scale,1,lag1,-1,dof,mother);

%%绘图      

figure(3);

%---绘制小波能量谱

subplot('position',[0.10.37 0.6 0.28])

levels =[0.0625,0.125,0.25,0.5] ;

Yticks =2.^(fix(log2(min(period))):fix(log2(max(period))));

contour(time,log2(period),log2(power),log2(levels));

xlabel('年份')

ylabel('周期/a')

title('a) 小波能量谱')

set(gca,'XLim',xlim(:))

set(gca,'YLim',log2([min(period),max(period)]),...

         'YDir','default', ...

         'YTick',log2(Yticks(:)), ...

         'YTickLabel',Yticks)

imagesc(time,log2(period),log2(power));

xlabel('年份')

ylabel('周期/a')

title('小波能量谱')

set(gca,'XLim',xlim(:))

set(gca,'YLim',log2([min(period),32]),...

         'YDir','default', ...

         'YTick',log2(Yticks(:)), ...

         'YTickLabel',Yticks)

%95% singificance contour,levelsat -99(fake)and 1(95% signif)

hold on

contour(time,log2(period),sig95,[-99,1],'k','linewidth',1);

hold on

plot(time,log2(coi),'k','linewidth',1)

hold off

%绘制全域能谱

subplot('position',[0.770.37 0.2 0.28])

plot(global_ws,log2(period))

hold on

plot(global_signif,log2(period),'--')

hold off

xlabel('能量谱')

title('全域能量谱')

set(gca,'YLim',log2([min(period),32]),...

         'YDir','default', ...

         'YTick',log2(Yticks(:)), ...

         'YTickLabel',Yticks)

set(gca,'XLim',[0,1.25*max(global_ws)])

画出来的图

dongyi

谢谢分享经验，总结的很好

丛薄

谢谢分享，继续学习一下

一大碗年糕

谢谢分享！

13top

很详细的总结了

苏拉苏拉

这么细致，棒棒的

棉花糖grace

楼主请问第3步中最大尺度为32是在哪里设置呀

vihenry2004

棉花糖grace 发表于 2021-7-6 16:51
楼主请问第3步中最大尺度为32是在哪里设置呀

Max.（<=32）这个选项

竹子

感谢楼主的经验分享 总结的很棒

Sibyl1208

感谢分享，最近一直在学习！很棒的帖子