MATLAB cwt连续小波变换 绘制小波振幅谱（影响锥coi）的三个函数示例（cwt contourf pcolor）

part 0 构造时间序列，并使用cwt函数计算小波系数矩阵

dt = 1;   % 一小时一个数
t = 0:dt:499;
Fs = 1/dt;
y = [10*sin(pi*(0:200)/5),5*cos(pi*(201:400)/8),1*sin(pi*(401:499)/10)];
y = y + 1*rand(size(t));
zq = 2:2:20;                                % 标注的周期

%% 以下amorlet小波基为例 计算小波系数
%[wt,f,coi] = cwt(___,fs) returns the cone of influence, coi, ... The cone of influence for the CWT is in Hz. 
%[wt,period,coi] = cwt(___,ts) returns the cone of influence, coi, ... The cone of influence for the CWT is in periods. 
[wt_y,f_y,coi_y] = cwt( y ,'amor',1/3600 );          %式一 FS = 1/3600 Hz → T = 3600s
[wt_Y,f_Y,coi_Y] = cwt( y ,'amor',1 );               %式二 FS = 1 cph     → T = 1 hour = 3600s
[wt_yT,period_yT,coi_yT] = cwt( y ,'amor',hours(1)); %式三 TS = 1 h → FS = 1 cph = 1/3600 Hz

T_y = 1./(3600*f_y);  % 式一输出的采样频率对应的周期（由赫兹Hz转换为小时对应周期）
T_Y = 1./f_Y;         % 式二输出的采样频率对应的周期（小时为单位）  
T_y-T_Y               % 尽管式一式二采样频率的单位不同，其输出频率对应的周期实质上应该一致
f_Y.*period_yT        % 两者应互为倒数，乘积为1小时，若式三为second(1),则乘积为1秒

%find(abs(wt_y-wt_Y)>1e-8)  % 式一与式三的小波系数实质上不一致（因为采样频率TS的数值不一致）有待解释

find(abs(wt_Y-wt_yT)>1e-8) % 式二与式三的小波系数实质上应该一致
% 以下方式二三以[式二计算结果 wt_Y f_Y coi_Y T_Y]绘制

diff(log2(T_Y))       % 60×1的列向量  并且所有值都是0.1

part 1 cwt 函数直接绘制

%% 绘制小波振幅谱的三种方式——之一      cwt不指定任何输出参数

figure()
cwt(y,'amor',hours(1));                    % cwt不指定任何输出参数，则直接绘图
set(gca,'YTickLabelMode','auto');
set(gca,'YTick',zq);
h = colorbar;
set(get(h,'Title'),'string','|\itW(a,b)|'); % colorbar标题
h.Label.String = 'magnitude';               % colorbar纵向标题
h.Label.String = '';                        % 清空纵向标题
title('Using cwt function plotting wavelet magnitude')

% pos = get(gcf,'position'); 
pos = [ 122         211        1655         449]; % 设置不同figure的大小，使其一致

part 2 contourf 绘制

%% 绘制小波振幅谱的三种方式——之二      使用contourf函数绘制

figure()
set(gcf,'position',pos)
contourf(t,log2(T_Y),abs(wt_Y),'linestyle','none')
colorbar
% set(gca,'YTickLabelMode','auto');
set(gca,'YTick',log2(zq));
set(gca,'YTicklabel',zq);
h = colorbar;
set(get(h,'Title'),'string','|\itW(a,b)|'); % colorbar标题
h.Label.String = 'magnitude';               % colorbar纵向标题
h.Label.String = '';                        % 清空纵向标题
title('Using contourf function plotting wavelet magnitude')
% plot coi
hold on
plot(t,-log2(coi_Y),'k--')
% plot(t,log2(1./coi_Y2),'k--') % 等价于上式

part 3 pcolor 绘制

%% 绘制小波振幅谱的三种方式——之三      使用pcolor函数绘制

figure()
set(gcf,'position',pos)
SP = pcolor(t,log2(T_Y),abs(wt_Y))
SP.EdgeColor = 'none';                      % 把不同矩阵位置的元素之间的间隔线设置为none
% SP.LineWidth = 0.0001;                    % pcolor属性只能是这样写，Egdecolor不对
SP.FaceColor = 'interp';                    % help pcolor 颜色更加连续
colorbar
% set(gca,'YTickLabelMode','auto');
set(gca,'YTick',log2(zq));
set(gca,'YTicklabel',zq);
title('Using pcolor function plotting wavelet magnitude')
% plot coi
hold on
plot(t,-log2(coi_Y),'k--')
% plot(t,log2(1./coi_Y),'k--') % 等价于上式

%% 注 方式二三所用coi的单位为频率（式二），标注为周期需取倒数，倒数关系在对数运算过程中变为负数关系，请参考cwt 的 文档， 以下是两种用法
%[wt,f,coi] = cwt(___,fs) returns the cone of influence, coi, ... The cone of influence for the CWT is in Hz. 
%[wt,period,coi] = cwt(___,ts) returns the cone of influence, coi, ... The cone of influence for the CWT is in periods. 

注意

这三种方式绘制的小波振幅谱黄色横线（振幅为10）与蓝色横线（振幅为5）的构造原始时间序列吻合较好

y = [10*sin(pi*(0:200)/5),5*cos(pi*(201:400)/8),1*sin(pi*(401:499)/10)];
y = y + 1*rand(size(t));
% 振幅为10 周期为10小时的振动在前201小时， 振幅为5 周期为16小时的振动在随后的200小时，振幅为1的正弦振动与随机噪声混叠，难以分辨

注意 cwt是直接画好小波振幅谱与coi， pcolor contourf 将coi之外的边缘部分，可将wt_Y的这部分设置为nan，则pcolor contourf的图这边缘部分将绘制为白色（不填充颜色）。