希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位画图

最后更新于：2023-11-16 22:10:55

一、代码运行环境

在使用该代码前，请务必安装时频域分析工具箱。（然后再在文末下载“希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位画图”相关代码并使用）

工具箱下载路径在这里

下载文件并解压，按照步骤操作即可完成工具箱安装。

已测试MATLAB2016/2018/2019，其他MATLAB版本理论上也可以使用。

推荐使用MATLAB2018a及更新版本。

二、文件说明

1.demoHHT.m

用于测试希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位功能的脚本文件，可以直接运行，其中调用了hhtSpec、marginalSpec、envSpec、InsFPA、pEMDandFFT。

2.hhtSpec.p

封装好的希尔伯特谱画图程序，p文件，可以调用。

函数参数说明：

function [hsfull,f,t] = hhtSpec(data,fs,type)
% 绘制信号希尔伯特谱
% 输入：
% data： 待分析信号
% fs：   采样频率，当fs未知时，可以将fs设置为空（[]），并将type设置为2。此时hht图纵轴将经过标准化
% type： 当type的值为1时，则优先采用MATLAB自带库中的hht函数进行画图，采用的绘图风格也与自带hht函数保持一致，此时hht函数的用法与MATLAB自带函数一致；
%        当type的值为2时，则强制采用第三方库文件中的希尔伯特谱函数进行画图，不过可能会存在内存不够无法顺利画图的可能
%        如果type没有传入参数，则该函数内将type设置为1
% 输出：
% hsfull：希尔伯特变换结果，即瞬时幅值
% f：频率轴数据
% t：时间轴数据

画图效果如下（参考）：

(1) 希尔伯特谱二维图

(2) 希尔伯特谱三维图（需使用完整版）

3.marginalSpec.p

封装好的边际谱计算和画图程序， p文件，可以调用。

函数参数说明：

function [mgS,f] = marginalSpec(imf,fs,type)
% 绘制边际谱
% 输入：
% imf：  imf分量，注意方向：imf是每行一个信号分量的矩阵
% fs：   采样频率
% type： 当type的值为1时，则优先采用MATLAB自带库中的函数进行运算；
%        当type的值为2时，则强制采用第三方库文件中的函数进行运算，不过可能会存在内存不够无法顺利画图的可能（数组超过预设的最大数组大小）
%        如果type没有传入参数，则该函数内将type设置为1
%        经测试，使用两种库函数做出的边际谱存在幅值差异，频率点基本一致
% 输出：
% mgS：  边际谱幅值。
% f：    边际谱的频率轴。
% 
% 典型用例：
% imf = kEMD(y);
% marginalSpec(imf,fs);

画图效果如下（参考）：

4.envSpec（完整版中获取）

封装好的包络谱计算和画图程序， p文件，可以调用。

函数参数说明：

function [envS,f,xEnv] = envSpec(data,fs)
% 求信号包络谱
% 输入：
% data： 待分析信号
% fs：   采样频率
% 输出：
% envS： 包络谱数值
% f：    包络谱频率轴
% xEnv： 包络线

画图效果如下（参考）：

5.InsFPA（完整版中获取）

封装好的瞬时频率/幅值/相位计算程序， p文件，可以调用。没有内置绘图程序，如需画图可以参考demoHHT.m文件中的写法。

函数参数说明：

function [insF,insP,insA] = InsFPA(data,fs)
% 求信号的瞬时频率、瞬时相位和瞬时幅值
% 输入：
% data： 目标信号
% fs：   采样频率
% 输出：
% insF： 瞬时频率
% insP： 瞬时相位
% insA： 瞬时幅值

6.kEMD.p

整合版EMD函数，整合了G-Rilling工具箱和MATLAB自带工具箱的EMD分解方法。

默认设置下按照MATLAB自带库进行。统一输出格式，配合后续希尔伯特-黄变换等运算无缝衔接。p文件，可以调用。

函数参数说明：

function imf = kEMD(data,type)
% 整合版EMD函数，整合了G-Rilling工具箱和MATLAB自带工具箱的EMD分解方法
% 默认设置下按照MATLAB自带库进行
% 并且统一输出格式，配合后续进行希尔伯特-黄变换等无缝衔接
% 输入：
% data：待分解的数据
% type： 如果type没有传入参数，例如调用时为kEMD(data)，则该函数内将type设置为1
%        当type的值为1时，则优先采用MATLAB自带库中的函数进行分解，此时emd函数的用法与MATLAB自带函数一致；
%        当type的值为2时，则强制采用第三方库文件中的函数进行分解
% 输出：
% 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)';

% 注意：在使用该代码之前，请务必安装工具箱：https://khsci.com/docs/index.php/2020/04/09/1/

7.pEMDandFFT.p

封装好的画图程序，可以绘制出信号EMD分解与各IMF分量频谱对照图。画图效果如下（参考）：

8.pEMD（完整版中获取）

封装好的EMD分解画图程序（不画频谱图），该程序在demoHHT中未用到。画图效果如下（参考）：

9.pFFT.p

封装好的fft计算程序，在pEMDandFFT.m中调用，一般不需要修改。

三、常见问题

暂无

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

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

如果需要封装好的希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位功能函数（hhtSpec.m、marginalSpec.m、envSpec.m、 InsFPA.m、kEMD.m、pEMDandFFT.m、pEMD.m和Fb_FFT.m）的无水印源码，可在下述连接获取。

此外如果想要得到包络谱和边际谱的计算结果数值也需要购买完整版代码哦~

编程不易，感谢支持~

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

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

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

获取通道一（淘宝）：点击此处获取完整版程序

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

（注意支付跳转失败的话，请使用浏览器打开本页面）

七、重要更新

20231115 （重要）增加了导出Hilbert变换结果（瞬时幅值、时间轴、频率轴）数据功能

20230305 （重要）增加了绘制Hilbert希尔伯特三维图的功能

20201220 修正了瞬时频率计算bug

20201009 完整版-修复了使用MATLAB2018时可能会报错的故障

20200819 代码初始版本（公开版所在版本）

八、常见问题

问题1：算例里面瞬时频率为何会有负值，该怎样处理？

在对 EMD 分解后的 IMF 分量求瞬时频率时，可能会出现负频率。这主要是因为在计算希尔伯特变换时，由于信号的局部特征或者数值计算误差，可能会导致相位计算的不稳定性。瞬时频率是相位的导数，因此这些不稳定性可能导致计算出的瞬时频率为负值。为解决这个问题，可以在计算瞬时频率时，对负值进行截断，即将负瞬时频率设为 0 或者用相邻的正瞬时频率值替换。

问题2：emd分解后的imf分量只能分别单独求瞬时频率嘛，能不能对完整的信号求瞬时频率？

许多实际信号都是非线性和非平稳的，EMD 可以较好地处理这类信号，EMD分解后的分量可以近似平稳，所以瞬时频率是对imf分量求的。对于如何对完整信号求整体瞬时频率，目前笔者未看到相关文献。

问题3：瞬时频率的单位是什么？

Hz（赫兹）。