最新EEMD介绍和使用

《最新EEMD介绍和使用》由会员分享，可在线阅读，更多相关《最新EEMD介绍和使用（22页珍藏版）》请在装配图网上搜索。

1、精品资料EEMD介绍和使用.EEMD介绍和使用1.什么是HHT？HHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。2.EMD分解的步骤。EMD分解的流程图如下：3.实例演示。给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10t)+5*sin(2*pi*35t)(1)为了对比，先用fft对求上述信号的幅频和相频曲线。1. function fftfenxi2. clear;clc;3. N=2048;4. %fft默认计

2、算的信号是从0开始的5. t=linspace(1,2,N);deta=t(2)-t(1);1/deta6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);7. % N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi;8. % x=(t=-200&t-200+N1*deta&t-200+N2*deta&t=200).*sin(w3*t);9. y = x;10. m=0:N-1;11. f=1./(N*deta)*m;%可以查看课本就是这样定义横坐标频率范围的12. %下面计算的Y就是x(t)的傅里叶变换数值13. %Y

3、=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移后-2,2之间的频谱值14. Y=fft(y);15. z=sqrt(Y.*conj(Y);16. plot(f(1:100),z(1:100);17. title(幅频曲线)18. xiangwei=angle(Y);19. figure(2)20. plot(f,xiangwei)21. title(相频曲线)22. figure(3)23. plot(t,y,r)24. %axis(-2,2,0,1.2)25. title(原始信号)复制代码（2）用Hilbert变换直接求该信号的瞬时频率

4、1. clear;clc;clf;2. %假设待分析的函数是z=t33. N=2048;4. %fft默认计算的信号是从0开始的5. t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);7. z=x;8. hx=hilbert(z);9. xr=real(hx);xi=imag(hx);10. %计算瞬时振幅11. sz=sqrt(xr.2+xi.2);12. %计算瞬时相位13. sx=angle(hx);14. %计算瞬时频率15. dt=diff(t);16. dx=diff(

5、sx);17. sp=dx./dt;18. plot(t(1:N-1),sp)19. title(瞬时频率)20.复制代码小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。Hilbert变换是求取瞬时频率的方法，但如果只用Hilbert变换求出来的瞬时频率也不准确。（出现负频，实际上负频没有意义！）（3）用HHT求取信号的时频谱与边际谱1. function HHT2. clear;clc;clf;3. N=2048;4. %fft默认计算的信号是从0开始的5. t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;6. x=5*sin(2*pi*1

6、0*t)+5*sin(2*pi*35*t);7. z=x;8. c=emd(z);9.10. %计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性11. m,n=size(c);12. for i=1:m;13. a=corrcoef(c(i,:),z);14. xg(i)=a(1,2);15. end16. xg;17.18. for i=1:m-119. %-20. %计算各IMF的方差贡献率21. %定义：方差为平方的均值减去均值的平方22. %均值的平方23. %imfp2=mean(c(i,:),2).224. %平方的均值25. %imf2p=mean(c(i,

7、:).2,2)26. %各个IMF的方差27. mse(i)=mean(c(i,:).2,2)-mean(c(i,:),2).2;28. end;29. mmse=sum(mse);30. for i=1:m-131. mse(i)=mean(c(i,:).2,2)-mean(c(i,:),2).2;32. %方差百分比，也就是方差贡献率33. mseb(i)=mse(i)/mmse*100;34. %显示各个IMF的方差和贡献率35. end;36. %画出每个IMF分量及最后一个剩余分量residual的图形37. figure(1)38. for i=1:m-139. disp(imf,

8、int2str(i) ;disp(mse(i) mseb(i);40. end;41. subplot(m+1,1,1)42. plot(t,z)43. set(gca,fontname,times New Roman)44. set(gca,fontsize,14.0)45. ylabel(signal,Amplitude)46.47. for i=1:m-148. subplot(m+1,1,i+1);49. set(gcf,color,w)50. plot(t,c(i,:),k)51. set(gca,fontname,times New Roman)52. set(gca,fontsi

9、ze,14.0)53. ylabel(imf,int2str(i)54. end55. subplot(m+1,1,m+1);56. set(gcf,color,w)57. plot(t,c(m,:),k)58. set(gca,fontname,times New Roman)59. set(gca,fontsize,14.0)60. ylabel(r,int2str(m-1)61.62. %画出每个IMF分量及剩余分量residual的幅频曲线63. figure(2)64. subplot(m+1,1,1)65. set(gcf,color,w)66. f,z=fftfenxi(t,z)

10、;67. plot(f,z,k)68. set(gca,fontname,times New Roman)69. set(gca,fontsize,14.0)70. ylabel(initial signal,int2str(m-1),Amplitude)71.72. for i=1:m-173. subplot(m+1,1,i+1);74. set(gcf,color,w)75. f,z=fftfenxi(t,c(i,:);76. plot(f,z,k)77. set(gca,fontname,times New Roman)78. set(gca,fontsize,14.0)79. yla

11、bel(imf,int2str(i),Amplitude)80. end81. subplot(m+1,1,m+1);82. set(gcf,color,w)83. f,z=fftfenxi(t,c(m,:);84. plot(f,z,k)85. set(gca,fontname,times New Roman)86. set(gca,fontsize,14.0)87. ylabel(r,int2str(m-1),Amplitude)88.89. hx=hilbert(z);90. xr=real(hx);xi=imag(hx);91. %计算瞬时振幅92. sz=sqrt(xr.2+xi.2

12、);93. %计算瞬时相位94. sx=angle(hx);95. %计算瞬时频率96. dt=diff(t);97. dx=diff(sx);98. sp=dx./dt;99. figure(6)100. plot(t(1:N-1),sp)101. title(瞬时频率)102.103. %计算HHT时频谱和边际谱104. A,fa,tt=hhspectrum(c);105. E,tt1=toimage(A,fa,tt,length(tt);106. figure(3)107. disp_hhs(E,tt1) %二维图显示HHT时频谱，E是求得的HHT谱108. pause109. figu

13、re(4)110. for i=1:size(c,1)111. faa=fa(i,:);112. FA,TT1=meshgrid(faa,tt1);%三维图显示HHT时频图113. surf(FA,TT1,E)114. title(HHT时频谱三维显示)115. hold on116. end117. hold off118. E=flipud(E);119. for k=1:size(E,1)120. bjp(k)=sum(E(k,:)*1/fs;121. end122. f=(1:N-2)/N*(fs/2);123. figure(5)124. plot(f,bjp);125. xlabe

14、l(频率 / Hz);126. ylabel(信号幅值);127. title(信号边际谱)%要求边际谱必须先对信号进行EMD分解128.129. function A,f,tt = hhspectrum(x,t,l,aff)130.131. error(nargchk(1,4,nargin);132.133. if nargin 2134.135. t=1:size(x,2);136.137. end138.139. if nargin 3140.141. l=1;142.143. end144.145. if nargin 4146.147. aff = 0;148.149. end150

15、.151. if min(size(x) = 1152. if size(x,2) = 1153. x = x;154. if nargin = 0214. error(inf doit etre 0)215. end216.217. M=max(max(im);218.219. im = log10(im/M+1e-300);220.221. inf=inf/10;222.223.224. imagesc(t,fliplr(1:size(im,1)/(2*size(im,1),im,inf,0);225. set(gca,YDir,normal)226. xlabel(time)227. y

16、label(normalized frequency)228. title(Hilbert-Huang spectrum)229. function f,z=fftfenxi(t,y)230. L=length(t);N=2nextpow2(L);231. %fft默认计算的信号是从0开始的232. t=linspace(t(1),t(L),N);deta=t(2)-t(1);233. m=0:N-1;234. f=1./(N*deta)*m;235. %下面计算的Y就是x(t)的傅里叶变换数值236. %Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi

17、*f)得到频移后-2,2之间的频谱值237. Y=fft(y);238. z=sqrt(Y.*conj(Y);239.240.241.242.复制代码4.总结。（1）边际谱与傅里叶谱的比较： 意义不同：边际谱从统计意义上表征了整组数据每个频率点的累积幅值分布，而傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分。 作用不同：边际谱可以处理非平稳信号，如果信号中存在某一频率的能量出现，就表示一定有该频率的振动波出现，也就是说，边际谱能比较准确地反映信号的实际频率成分。而傅里叶变换只能处理平稳信号。（2） HHT与Hilbert变换的比较： Hilbert变换只是单纯地求

18、信号的瞬时振幅，频率和相位，有可能出现没有意义的负频率；HHT变换先将信号进行EMD分解，得到的是各个不同尺度的分量，对每一个分量进行Hilbert变换后得到的是有实际意义的瞬时频率。PS:运行上面的程序需要装时频工具箱，我仅将用到的emd分解的程序贴到下面。把标准的eemd函数加载到set path 或者直接放在 运行程序的文件夹内之后，即可调用eemd函数。具体的函数使用如下。该函数为function allmode=eemd(Y,Nstd,NE)输入的英文注释如下：% INPUT:% Y: Inputted data;% Nstd: ratio of the standard devia

19、tion of the added noise and that of Y;% NE: Ensemble number for the EEMD即：Y是输入的数据；Nstd是用来设置添加高斯白噪声的标准差的，用以消去原信号中的噪声，Nstd要根据原信号中的噪声干扰大小具体情况而定，高斯白噪声的标准差设置一般为0.010.4,具体设置没有一个确定的公式，根据信号来确定；NE是用来设置添加噪声的次数，NE通常取50或100。输入的英文注释如下：% OUTPUT:% A matrix of N*(m+1) matrix, where N is the length of the input%data

20、 Y, and m=fix(log2(N)-1. Column 1 is the original data, columns 2, 3, .% m are the IMFs from high to low frequency, and comlumn (m+1) is the% residual (over all trend).即输出一个N*(M+1)的矩阵，其中N是输入数据Y的长度，m=fix(log2(N)-1。第一列是原始数据，之后的是IMFs。信号处理中，频率是信号最重要的表示。传统的傅里叶变换分析方法并不能分析出信号的某一频率在甚么时刻出现，为此产生了能同时在时间和频率上表示信

21、号密度和强度的时频分析，如短时傅里叶变换和小波变换等，但其基本思想都是根据傅里叶分析理论，对非线性非平稳信号的分析能力不足，受限于Heisenberg不确定原理。HHT ( Hilbert Huang Transform)是由N. E.Huang 等人在1998 年提出的一种崭新的时频分析方法，能够对非线性非平稳的信号进行分析，同时具有良好自适应性的特点。其本质是对信号进行平稳化处理，将具有不同时间尺度的信号逐级分解开来。HHT 方法在各领域已得到了广泛应用，但依然存在一些不足，例如易产生虚假分量和模态混叠等。针对传统经验模式( Empirical Mode Decomposit ion，EM

22、D)分解方法所导致的模态混叠现象，法国以Flandrin 为首的EMD 算法研究小组和Huang 本人的研究小组通过对EMD 分解白噪声结果统计特性的大量研究，提出通过加噪声辅助分析( NA DA ) 的EEMD ( EnsembleEmpirical Mode Decomposition) 方法，将白噪声加入信号来补充一些缺失的尺度，在信号分解中具有良好的表现。EEMD仿真系统的实现利用了Matlab 平台，通过GUI 控件实现了系统设计，能直观方便地进行比较分析，验证了EEMD 在抗混叠方面较原有方法的改进。1 经验模式分解( EMD) 和IMFHHT 方法包含两个主要步骤:( 1) 对原

23、始数据进行经验模式分解( EMD) ，把数据分解为满足Hilbert 变换要求的n 阶本征模式函数( IMF) 和残余函数之和。( 2) 对每一阶IMF 进行Hilbert 变换，得到瞬时频率，从而求得时频图。函数必须关于时间轴局部对称，且其过零点与极值点个数相同。此类函数被称为固有模态函数( Int rinsicMode Function，IMF) 。经验模式分解方法能把非平稳、非线性信号分解成一组稳态和线性的序列集，即本征模式函数。根据Huang 的定义，每一阶的IMF 应满足两个条件:( 1) 数据的极值点和过零点交替出现，且数目相等或最多相差一个任何点上；( 2) 在任何点上，有局部最

24、大值和局部最小值定义的包络的均值必须是零。其筛选算法如下:( 1) 对于输入信号x ( t) ，确定x ( t) 所有极值点。( 2) 用三次样条函数对极大点和极小点分别进行拟合得到x ( t) 的上下包络线。( 3) 用原始数据序列减去上下包络线的均值。平均曲线：细节信号：( 4) 通常s( t ) 还不满足IMF 的条件，需重复进行以上步骤，进行迭代处理，H uang 给出的迭代停止准则为：SD 是筛选门限值，一般取值为0. 2 0. 3，若计算SD 小于这个门限值，筛选迭代将会结束。经过n 次迭代满足停止准则后得到的sn ( t) 即为有效IMF，剩余信号则进入下一轮筛选过程。经过多次筛

25、选后，原始数据序列被分解为一组IMF 分量和一个残余量，得到的IMF 都是平稳的，通过Hilbert 变换得到的结果能够很好地分析非线性非平稳的信号。2 传统EMD 的不足与缺陷当信号的时间尺度存在跳跃性变化时，对信号进行EMD 分解，会出现一个IMF 分量包含不同时间尺度特征成分的情况，称之为模态混叠。模态混叠的出现一方面和EMD 的算法有关，另一方面也受原始信号频率特征的影响。Huang 曾经提出了中断检测的方法来解决模态混叠现象，即直接对结果进行观察，如果出现混叠则重新分解，这种方法需要人为后验判断。重庆大学的谭善文提出了多的EMD 思想，对每一个IMF 规定一个尺度范围来解决模态混叠，

26、但是这种方法牺牲了EMD 良好的自适应性。3 引入正态分布白噪声的EEMD为了更好地解决模态混叠问题，Huang 提出了EEMD，这是一种噪声辅助信号处理方法。降噪技术的目的是将噪声从信号中去除，不过在一些情况下，可以通过加入噪声的方法来进行辅助分析,这钟方法就称为噪声辅助信号处理( NADA) ，噪声辅助信号处理方法最常见的就是预白化。在信号中加入白噪声来平滑脉冲干扰，被广泛用于各种信号分析领域。在EMD 方法中，得到合理IMF 的能力取决于信号极值点的分布情况，如果信号极值点分布不均匀，会出现模态混叠的情况。为此，Huang 将白噪声加入待分解信号，利用白噪声频谱的均匀分布，当信号加在遍布

27、整个时频空间分布一致的白噪声背景上时，不同时间尺度的信号会自动分布到合适的参考尺度上，并且由于零均值噪声的特性，经过多次平均后，噪声将相互抵消,集成均值的结果就可作为最终结果。EEMD 步骤如下:( 1) 向信号加入正态分布白噪声。( 2) 将加入白噪声的信号分解成各IMF 分量。( 3) 重复步骤( 1) ，( 2) ，每次加入新的白噪声序列。( 4) 将每次得到的IMF 集成均值作为最终结果。EMMD 算法流程如图1 所示。图1 EEMD 算法流程图4 系统功能介绍和仿真实验分析为了验证EEMD 方法的改进之处，利用Mat lab 的GU I 工具设计了简单直观的仿真系统。此系统实现的功能

28、是，对输入信号进行传统EMD分解和EEMD 分解，可显示信号分解后的各个模态函数IMF 分量及其瞬时频率，并能对Hilbert 时频谱进行刻画。系统界面如图2 所示。图2 仿真系统界面参数设置功能 可自由设置加入白噪声的方差和噪声组数目( 范围1 500) ，当方差设置为0，噪声组数目选择为1 时，该系统实现传统EMD 分解的功能。EEMD 分解功能 对信号进行加入上述设定白噪声EEMD 分解，并刻画出输入信号的Hilbert 时频谱。显示IMFs 功能 可通过弹出FIG 的形式显示对信号分解后的各IMF 分量及瞬时频率。仿真实验结果如下：首先对多分量理想样本信号进行分解，信号构成如下：其中，

29、归一化频率为:EMD 分解方法应将包含4 个频率分量的信号分解为4 个包含单一频率信息的IMF 分量。分解结果如图3 所示。图3 传统EMD 对理想信号H ilber t 谱图可以看到，对于无干扰的理想信号，传统EMD 分解方法具有非常好的效果，清晰地将4 个频率分量在Hilbert 谱上显示了出来。对一组存在中断干扰的实际信号进行分解，结果如图4 图6 所示。图4 实际信号时域图图5 传统EMD 对信号的分解图6 传统EM D 对信号的H ilber t 谱刻画通过频谱图可以看到，低频分量混杂在一起，难以分辨。对EEMD 分解方法进行分析，加入了100 组标准差为0. 2 的高斯白噪声，结果如图7，图8 所示。通过Hilbert 谱的比较可以看出，分解结果有了较大改进。图7 EEMD 对信号的分解图8 EEMD 对信号的H ilber t 谱刻画5 结 语EEMD 以噪声辅助信号处理原理为基础，通过加入小幅度的白噪声来均衡信号，有效地解决了模态混叠现象，利用高斯白噪声零均值的特性，使真实信号得到了保留，是对传统EMD 分析方法的巨大改进。