MATLAB自带函数实现经验模态分解
MATLAB从2018a开始给出了内置函数来实现经验模式分解(EMD)与希尔伯特-黄变换(HHT),函数名分别是emd与hht,用户可以直接调用两个函数来实现经验模式分解以及与希尔伯特-黄变换。
Syntax(语法)
-
[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和残差信号相加就是你输入的那个信号。
-
[imf,residual,info] = emd(X)
returns additional information info on IMFs and residual signal for diagnostic purposes.
-
[] = emd(,Name,Value)
performs the empirical mode decomposition with additional options specified by one or more Name,Value pair arguments. 类似键值对的输入方式,先输入键名,再输入值。
-
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选择器。