MATLAB自带函数实现经验模态分解总结

MATLAB自带函数实现经验模态分解

MATLAB从2018a开始给出了内置函数来实现经验模式分解（EMD）与希尔伯特-黄变换（HHT），函数名分别是emd与hht，用户可以直接调用两个函数来实现经验模式分解以及与希尔伯特-黄变换。

Syntax(语法)

1. [imf,residual] = emd(X)

    returns intrinsic mode functions imf and residual signal residual corresponding to the empirical mode decomposition of X. Use emd to decompose and simplify complicated signals into a finite number of intrinsic mode functions required to perform Hilbert spectral analysis.
   
    意思就是返回值是一个经过emd处理后的本征函数imf和残差信号。验证一下，把imf和残差信号相加就是你输入的那个信号。
   

2. [imf,residual,info] = emd(X)

    returns additional information info on IMFs and residual signal for diagnostic purposes.
   

3. [] = emd(,Name,Value)

    performs the empirical mode decomposition with additional options specified by one or more Name,Value pair arguments.
    	
    类似键值对的输入方式，先输入键名，再输入值。
   

4. emd(___)

    plots the original signal, IMFs, and residual signal as subplots in the same figure.
   

Matlab里面的例子

1.对信号进行经验模态分解和希尔伯特谱可视化
load(‘sinusoidalSignalExampleData.mat’,‘X’,‘fs’)
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel(‘Time(s)’)

sinusoidalSignalExampleData.mat这个数据在matlab文件里，使用命令help emd，然后进入说明文档那里，点开OpenLive Script就有了。

混合信号包含不同振幅和频率值的正弦波。

要创建希尔伯特频谱图，需要信号的固有模态函数(IMFs)。进行经验模态分解，计算信号的imf和残差。由于信号不是平滑的，指定’pchip’作为插值方法。

[imf,residual,info] = emd(X,‘Interpolation’,‘pchip’);

使用上面那条命令的时候，可能会报错：bad option field name: Interpolation，原因是因为你安装了EMD工具箱，这与matlab的emd函数重复导致的。只要把EMD工具箱先移出matlab的路径就可以了。

这条命令出来的结果是多了一个info的1×1struct，里面有6个字段，分别是NumIMF，NumExtrema,NumZerocrossing,NumSifting,MeanEnvelopeEnergy和RelativeTolerance。

使用经验模态分解得到的imf分量创建希尔伯特谱图。 命令：hht(imf,fs)

频率与时间的关系图是一个稀疏的图，其中有一个竖直的彩色条，表示IMF中每个点的瞬时能量。图中是原始混合信号分解后各分量的瞬时频谱。三个imf出现在图中，其频率在1秒内发生明显变化。

2.正弦波本征模态函数的过零和极值

生成两个正弦波s和z，这样s是三个正弦波的和，z是一个调幅的单一正弦波。通过计算两个信号差的无穷范数来验证两个信号是相等的。

t = 0:1e-3:10;
omega1 = 2*pi*100;
omega2 = 2*pi*20;
s = 0.25*cos((omega1-omega2)*t) + 2.5*cos(omega1*t) + 0.25*cos((omega1+omega2)*t);
z = (2+cos(omega2/2*t).^2).*cos(omega1*t);
norm(s-z,Inf) 

ans = 3.2729e-13

绘制正弦曲线并选择从2秒开始的1秒间隔。

plot(t,[s' z'])
xlim([2 3]) 
xlabel('Time (s)')
ylabel('Signal')

获取信号的声谱图。谱图显示了三个不同的正弦分量。傅里叶分析把这些信号看作是正弦波的叠加。

使用emd来计算信号的固有模态函数(IMFs)和附加的诊断信息。默认情况下，该函数输出一个表，其中指示每个IMF的筛选迭代次数、相对容忍度和筛选停止条件。经验模态分解把信号看成z。

零交叉和局部极值的数目最多相差1。这满足了信号为IMF的必要条件。

绘制IMF曲线，并选择从2秒开始的0.5秒间隔。IMF是一个AM信号，因为emd认为信号是调幅的

plot(t,imf)
xlim([2 2.5])
xlabel('Time (s)')
ylabel('IMF')

3.计算振动信号的固有模态函数
略（目前好像用不着）

4.可视化信号的残差和本征模态函数
加载并可视化一个由频率变化明显的正弦波组成的非平稳连续信号。电钻的振动和烟花的声音都是非平稳连续信号的例子。

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

混合信号包含不同振幅和频率值的正弦波。

进行经验模态分解，绘制信号的本征模态函数和残差。由于信号不是平滑的，指定’pchip’作为插值方法。

命令：emd(X,‘Interpolation’,‘pchip’,‘Display’,1)

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
    1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
    2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
    3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
    4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
    5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
    6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
    7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
    8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
    9      |        2     |      0.15923   |  SiftMaxRelativeTolerance

The decomposition stopped because the number of extrema of the residual signal is less than ‘MaxNumExtrema’.

上面这个命令给出了一个迭代次数的表格，还有一个Figure：经验模态分解的各imf分量，但是图中只显示3个imf分量和1个残差，如果要查看所有imf分量，就要用鼠标在途中右键，选择imf选择器。