IMF方差贡献率、平均周期、相关系数的计算,高频、低频、趋势项分量判别与重构程序(公开版)

最后更新于:2023-06-18 15:42:42

一、运行环境

MATLAB2018a或更新版本已经过测试,之前的版本未做测试。

*注:本程序中使用EEMD作为测试方法,但imfHLdif和imfClc程序对EMD、CEEMD、CEEMDAN等“类EMD”方法都可以适用(VMD和EWT不适用于高、低频判别和重构),但是就需要同学们自行在程序中做适应性替换啦。

二、文件说明

注:图标代表该m文件为脚本文件,可以直接运行;图标代表函数文件,在没有输入变量的情况下无法直接运行。更详细的解释可以看这里

1.testIMF.m

用于测试计算IMF方差贡献率、平均周期、相关系数,高频、低频、趋势项分量判别与重构功能的脚本文件,可以直接运行,其中调用了imfClc和imfHLdif,演示了这两个函数的调用方法,以及在命令行窗口打印出计算结果;原始数据图;EEMD分解结果图;绘制高频、低频、趋势项图的用法。运行后会得到下述结果:

方差贡献率、平均周期、相关系数、高频和低频判断等的计算结果

2. imfClc.p

封装好的方差贡献率、平均周期、相关系数的计算程序,p文件。

函数参数说明:

function [VarR,AvePer,PearsonCor] = imfClc(data,imf)
%%
% 对于“类EMD”方法分解后得到的各个分量计算评价指标
% 包括方差贡献率、平均周期和Pearson相关系数
% 输入:
% data:分解前的原数据
% imf:分解后得到的分量,注意imf需要沿行方向分布
% 输出:
% VarR:方差贡献率
% AvePer:平均周期
% PearsonCor:Pearson相关系数

3.imfHLdif.p

封装好的程序,可以根据重构算法将分解得出的IMF进行高低频的区分,内置了绘制 高频、低频、趋势项分解图的方法,不过调用该函数画图,图片的横坐标为数据点数,而不是日期,如果要以日期为横坐标画图,需要参考 testIMF 脚本的写法。p文件。

函数参数说明:

function [HighCom,LowCom,TrCom,HighIdx,LowIdx]=imfHLdif(data,imf,figflag)
%% 根据重构算法将分解得出的IMF进行高低频的区分
% 参考《基于EEMD模型的中国碳市场价格形成机制研究》
% 该方法将IMF1记为指标1,IMF1+IMF2为指标2,以此类推,
% 前i个IMF的和加成为指标i,并对该均值是否显著区别于0进行t检验。
% 输入:
% data:分解前的原始数据
% imf:经过模态分解方法得到的分量,每一行为一个分量
% figflag:设置是否画图的参数,'on'为画图,'off'为不画图
% 输出:
% HighCom:重构后的高频分量
% LowCom:重构后的低频分量
% TrCom:趋势项
% HighIdx:高频分量对索引
% LowIdx:低频分量对索引

画图效果如下(参考):

4.pEEMD.p

封装好的EEMD画图程序,其中调用了kEEMD.p,p文件,可以调用。

函数参数说明:

function  imf = pEEMD(data,FsOrT,Nstd,NE)
% 画信号EEMD分解图
% 输入:
% y为待分解信号
% FsOrT为采样频率或采样时间向量,如果为采样频率,该变量输入单个值;如果为时间向量,该变量为与y相同长度的一维向量。如果未知采样频率,可设置为1
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% 输出:
% imf为经eemd分解后的各imf分量值
% 例1:(FsOrT为采样频率)
% fs = 100;
% t = 1/fs:1/fs:1;
% data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pEEMD(data,fs,0.2,100);
% 例2:(FsOrT为时间向量,需要注意此时FsOrT的长度要与y相同)
% t = 0:0.01:1;
% data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pEEMD(data,t,0.2,100);

5.kEEMD.p

调整后的EEMD函数,整合了中央大学工具箱的EEMD分解方法。

调整包括:

1.将imf分量由纵向改为横向,使之与时频域分析工具箱中其他模态分解工具的格式保持一致。方便与后续方法的衔接

2.原始eemd函数中第一个分量为原始信号,删之

函数参数说明:

function imf = kEEMD(data,Nstd,NE)
% 调整后的EEMD函数,整合了中央大学工具箱的EEMD分解方法
% 调整包括:1.将imf分量由纵向改为横向,使之与时频域分析工具箱(http://hk.khscience.cn/docs/index.php/2020/04/09/1/)中其他模态分解工具的格式保持一致。方便与后续方法的衔接
%          2.原始eemd函数中第一个分量为原始信号,删之
% 输入:
% data:待分解的数据
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% 输出:
% imf:内涵模态分量,统一为n*m格式,其中n为模态数,m为数据点数。例如 imf(1,:)即IMF1,imf(emd,:)即为残差

% 注意:使用两种库所得到的imf结果会存在差异,这是由两个库自身算法差异导致的,属于正常现象。
% 有余力的同学可以根据自己的需要修改下述代码
% 例如当type设置为1时,mat_emd的输入参数可以完全参照下述文档进行设置:https://ww2.mathworks.cn/help/signal/ref/emd.html?s_tid=srchtitle
% 比如设置最大IMF数为5,则①行的代码可以写为:imf = mat_emd(data,'MaxNumIMF',5)';

% 注意:在使用该代码之前,请务必安装工具箱:http://hk.khscience.cn/docs/index.php/2020/04/09/1/

6.eemd.m和extrema.m

第三方EMD工具箱里的开源程序,为了方便使用整合到一起了。这两个代码为开源,不属于出售和售后范围,它们也可以在时频域工具箱中免费获取到。

三、使用教程

1.运行测试脚本

先在MATLAB里打开下载好的文件夹,然后运行 testIMF.m程序,程序运行完毕后如果没报错,则说明运行环境正常,程序正确。

2.导入数据

复制一个testIMF.m的文件副本,在这个副本里做如下修改:

根据你的文件类型的不同(excel,txt,csv等),将数据导入MATLAB的方法有所不同。同学们可以先参考这个文档。或者看博主针对常用文件的导入方法的这个教程

如果你已经实现了数据导入,这时候应该拥有了一个一维数据变量,这时候就可以调用函数进行计算了。

3.实现特征提取

参照 testIMF.m文件中1-5步, imfClcimfHLdif函数的调用方法以及在命令行打印结果值、绘图方法,运行程序即可。

四、 关于完整版与公开版代码

公开版代码的函数文件为p文件,可以被调用,但无法查看代码。完整版代码中全部为m文件,m文件可以查看源码并自由修改。

公开版代码仅可对长度300以内的数据分解;完整版无长度限制。

完整版代码画图无水印。

如果需要封装好的画图函数( testIMF.m、 imfClc.m、 imfHLdif.m、 pEEMD.mkEEMD.m)的源码,可在下述连接(完整版)获取。

编程不易,感谢支持~

五、获取公开版程序(需使用电脑浏览器打开)

IMF数据指标计算与高低频判断(公开版)

注:公开版代码需使用MATLAB2022a及以上版本

六、获取完整版程序(使用电脑浏览器或者手机浏览器打开)

获取通道一(淘宝):点击此处获取完整版程序

获取通道二(本页面):点击下面“立即支付”按钮,付款后获取完整版代码下载链接和售后联系方式~本通道处于测试阶段,使用该通道可以额外优惠(仅需38元)。付款完成后刷新一下本页面即可看到下载链接。

(注意支付跳转失败的话,请使用浏览器打开本页面。)(付款完成后刷新本页面就能看到下载链接)

您需要先支付 38元 才能查看此处内容!立即支付

七、完整版代码重要更新

20230618 绘制图片横坐标由日期改为采样点,更具有普适性。

20220520 修复了当单个imf分量被判定为高频分量或低频分量时程序会报错的bug。

20220320 代码初始版本