EMD分解的流程图如下
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默认计算的信号是从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).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(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=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变换直接求该信号的瞬时频率1.clear;clc;clf;2.%假设待分析的函数是z=t^33.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(sx);17.sp=dx./dt;18.plot(t(1:N-1),sp)19.title('瞬时频率')20.复制代码小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。
EMD经验模式分解信息汇总资料
EMDEmpirical Mode Decomposition经验模态分解美国工程院院士黄锷1998年提出一种自适应数据处理或挖掘方法，适用于非线性、非平稳时间序列的处理。
1.什么是平稳和非平稳时间序列的平稳，一般是宽平稳，即时间序列的方差和均值是和时间无关的常数，协方差与与时间间隔有关、与时间无关。
未来样本时间序列，其均值、方差、协方差必定与已经获得的样本相同，理解为平稳的时间序列是有规律且可预测的，样本拟合曲线的形态具有“惯性”。
而非平稳信号样本的本质特征只存在于信号所发生的当下，不会延续到未来，不可预测。
严格来说实际上不存在理想平稳序列，实际情况下都是非平稳。
2.什么是EMD经验模态分解方法？EMD理论上可以应用于任何类型时间序列信号的分解，在实际工况中大量非平稳信号数据的处理上具有明显优势。
这种优势是相对于建立在先验性假设的谐波基函数上的傅里叶分解和小波基函数上的小波分解而言的。
EMD分解信号不需要事先预定或强制给定基函数，而是依赖信号本身特征自适应地进行分解。
相对于小波分解：EMD克服了基函数无自适应性的问题，小波分析需要选定一个已经定义好的小波基，小波基的选择至关重要，一旦选定，在整个分析过程中无法更换。
这就导致全局最优的小波基在局部的表现可能并不好，缺乏适应性。
而EMD不需要做预先的分析与研究，可以直接开始分解，不需要人为的设置和干预。
相对于傅里叶变换：EMD克服了传统傅里叶变换中用无意义的谐波分量来表示非线性、非平稳信号的缺点，并且可以得到极高的时频分辨率。
EMD方法的关键是将复杂信号分解为有限个本征模函数IMF，Intrinsic Mode Function。
分解出来的IMF分量包含了原信号的不同时间尺度上的局部特征信号。
这句话中：不同时间尺度=局部平稳化，通过数据的特征时间尺度来获得本征波动模式，然后分解or筛选数据。
本质上，EMD将一个频率不规则的波化为多个单一频率的波+残波的形式。
希尔伯特-黄变换在地质雷达资料处理中的应用
Hilbert Huang Transform and Applications of it in Ground Penetrating Radar Data Analysis and Processing
ABSTRACT
Hilbert-Huang Transform is an algorithm which apply huang transform and hilbert transform to original signal in proper order. It applys in Ground Penetrating Radar data analysis and processing due to its time frequency analysis advantages in order that (1)we can depict the Ground Penetrating Radar record wave characteristics much more finely in time frequency domain to obtain particular and benefit strata information;(2) we can get rid of noise in noise contained Ground Penetrating Radar data to obtain better quanlity,according to the decomposition advantages of huang transform. (3) we can extract natural Ground Penetrating Radar attribute profiles to verify or distinguish the wave impedance and tiny structures,according to the Ground Penetrating Radar data and the work flow for instantaneous attributes. The research method in this thesis is Hilbert-Huang Transform which is a combination of huang transform and hilbert transform.Huang transform can decompose an original signal to a series of intrinsic mode functions which its number is limited.We can get the meaningful instantaneous frequency to display exquisitely in time frequency domain after applying hilbert transform to IMFs. Many numerical algorithms correspond to Hilbert Huang Transform are stimulated and applied in actual Ground Penetrating Radar data processing due to its advantages and characteristics.Huang transform is a superior method because components decomposed by huang transform is much more pure than that of wavelet transform decomposed.The“end effect” phenomena are better suppressed after improving the
EMD分解的流程图如下
1. 什么是HHTHHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。
2. EMD分解的步骤EMD分解的流程图如下:3. 实例演示。
给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5si n(2*pi*10t)+5*s in (2*pi*35t)(1)为了对比，先用fft对求上述信号的幅频和相频曲线。
1. function fftfenxi2. clear;clc;3. N=2048;4. %fft默认计算的信号是从0开始的5. t=li nspace(1,2,N);deta=t(2)-t(1);1/deta6. x=5*si n(2*pi*10*t)+5*si n( 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).*si n(w1*t)+(t>-200+N1*deta &t<=-200+N2*deta).*si n(w2*t)+(t>-200+N2*deta&t<=200).*si n(w3*t);9. y = x;10. m=0:N-1;11. f=1./(N*deta)*m;% 可以查看课本就是这样定义横坐标频率范围的12. 虾面计算的Y就是x(t)的傅里叶变换数值13. %Y=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变换直接求该信号的瞬时频率1. clear;clc;clf;2. 液设待分析的函数是z=t A33. N=2048;4. %fft默认计算的信号是从0开始的5. t=li nspace(1,2,N);deta=t(2)-t(1);fs=1/deta;6. x=5*si n(2*pi*10*t)+5*si n( 2*pi*35*t);7. z=x;8. hx=hilbert (z);9. xr=real(hx);xi=imag(hx);10. %十算瞬时振幅11. sz=sqrt(xr.A2+xi.A2);12. %十算瞬时相位13. sx=angle(hx);14. %十算瞬时频率15. dt=diff(t);16. dx=diff(sx);17. sp=dx./dt;18. plot(t(1:N-1),sp)19. title(' 瞬时频率')20.复制代码小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。
EMD详解及其应用
EMD 详解及其应用王骏一G2*******EMD ，全称Eeath Movers' Distance ，它是用来衡量两个特征分布之间相似度的一个重要的度量，在我们的科研工作中起到了相当重要的作用。
EMD 的前身——运输问题运输模型是指，设某种物资有m 个产地A 1,A 2，…,A m ，供应量分别为a 1，a 2，…，a m 个单位；联合供应n 个销地B 1，B 2,…,B n ,需求量分别为b 1，b 2,…,b n 个单位（总供应量大于等于总需求量）。
假设从产地A i 向销地B j 运输一个单位物资的费用为c ij ，怎样调运物资才能使运输费用最少？记从产地A i 到销地B j 的运输量为x ij ，则总的运输成本可记为：∑===n 1j 1i ijij x ×c S ，我们的目标是求出S 的最小值，即min （S ）。
运输问题表格EMD 借用了运输问题求解的思路，它可以被理解为“从一种分布变换为另一种分布的最小代价”，它最早被Peleg ，Werman 和Rom 介绍应用于计算视觉问题。
后来，人们将该流程移植到特征分布的比较中，把一个特征分布当作“供货商”，而另一个为“消费商”。
2ij 1ij P C P C 特征分布特征分布消费商供货商−→−−→−定义C ij 为从第一个特征分布的第i 个元素与第二个个特征分布的第j 个元素之间的“距离”（C ij 可以是任何距离的度量，应根据当前处理的问题灵活选择）。
再使用运输问题的算法找到最优路径矩阵，就得到两个特征分布之间的EMD 。
设两个区域RA 和RB ，可用区域内某一特征信息的概率分布分别表征为:{(rA1,vA1),(rA2,vA2),...,(rAm,vAm )}{(rB1,vB1),(rB2,vB2),...,(rBn ,vBn )}则区域RA 和RB 的EMD 可以定义为:货物概率分布直方图消费商特征分布供货商特征分布→→→21P P ⎪⎪⎪⎩⎪⎪⎪⎨⎧∑=∑∑≤≤≤≤≤≤≤≤≤≤≥∑∑∑∑=∑∑∑==========)min(),(1,),(1,),(1,1,0),(f S.T.min )R ,EMD(R 1,m 1i n 111),(),(),(B A 1m 1i n 1m 1i n 1m i Bj Ai j m i Bj n j Ai j i f j i d j i f n j j j v v j i f n j v j i f m i v j i f n j m i j i我们以这两个集合为例：{'a','a','a','b','b','c','d','d','d','d'}和{'a','a','c','c','c','c','c','e','k','k'}，他们的概率密度分布图依次为：根据概率密度分布图，我们可以得到EMD 转化表格（如下图）：EMD(str1,str2)=2*0+2*1+1*0+2*1+1*1+1*10+1*7=230.00%10.00%20.00%30.00%40.00%50.00%a b c d 0%10%20%30%40%50%a c e k我们称EMD(str1,str2)=23为基本可行解，进而，我们利用“表上作业法”求出最优解：EMDmin(str1,str2)。
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默认计算的信号是从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).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(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=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变换直接求该信号的瞬时频率1.clear;clc;clf;2.%假设待分析的函数是z=t^33.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(sx);17.sp=dx./dt;18.plot(t(1:N-1),sp)19.title('瞬时频率')20.复制代码小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。
用希尔伯特黄变换(HHT)求时频谱和边际谱
用希尔伯特黄变换（HHT）求时频谱和边际谱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对求上述信号的幅频和相频曲线。
代码:function fftfenxiclear;clc;N=2048;%fft默认计算的信号是从0开始的t=linspace(1,2,N);deta=t(2)-t(1);1/detax=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);% N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi;%x=(t>=-200&t<=-200+N1*deta).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*det a).*sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t);y = x;m=0:N-1;f=1./(N*deta)*m;%可以查看课本就是这样定义横坐标频率范围的%下面计算的Y就是x(t)的傅里叶变换数值%Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移后[-2,2]之间的频谱值Y=fft(y);z=sqrt(Y.*conj(Y));plot(f(1:100),z(1:100));title('幅频曲线')xiangwei=angle(Y);figure(2)plot(f,xiangwei)title('相频曲线')figure(3)plot(t,y,'r')%axis([-2,2,0,1.2])title('原始信号')（2）用Hilbert 变换直接求该信号的瞬时频率代码:clear;clc;clf;%假设待分析的函数是z=t^3N=2048;%fft 默认计算的信号是从0开始的t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);z=x;hx=hilbert(z);xr=real(hx);xi=imag(hx);%计算瞬时振幅sz=sqrt(xr.^2+xi.^2);%计算瞬时相位sx=angle(hx);%计算瞬时频率dt=diff(t);dx=diff(sx);sp=dx./dt;plot(t(1:N-1),sp)title('瞬时频率')小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。
经验模态分解算法
经验模态分解算法
EMD算法的步骤如下：
1.将要分解的信号称为原始信号，记为x(t)。
2.寻找x(t)的极大值点和极小值点，这些点将原始信号分为一系列小段。
3.对每个小段进行插值，使均匀分布的数据点可以拟合出这个小段。
4. 利用Cubic Spline插值法或其他插值方法找到一个包络线，该包络线连接这些插值点的极大值点和极小值点。
即为信号中的一条上包络线和一条下包络线。
5.计算出平均值函数m(t)=(上包络线+下包络线)/2
6.计算x(t)与m(t)的差值d(t)=x(t)-m(t)。
7.如果d(t)是一条IMF，则终止算法；否则将d(t)作为新的原始信号，重复步骤2-6
8.将计算出的IMF组合起来，得到原始信号x(t)的EMD分解结果。
EMD算法的特点是对信号进行自适应分解，能够捕捉到不同频率的局部特征。
它不需要提前设定基函数或者滤波器，而是根据信号中的局部特征自动适应地生成各个IMF。
因此，EMD算法在信号处理领域中得到了广泛应用，如地震信号分析、生物信号处理等。
然而，EMD算法也存在一些问题。
其中最主要的问题是固有模态函数的提取过程中可能出现模态混叠的情况，即两个或多个IMF的频率相似且在一些区间内相互重叠，使得提取的IMF不纯粹。
为了克服这个问题，研
究者们提出了一些改进的EMD算法，如快速EMD、改进的EMD等。
这些改进方法在一定程度上提高了EMD算法的可靠性和稳定性。
总之，经验模态分解算法是一种有效的信号分解方法，能够提供信号的局部特征表示。
它在很多领域有广泛的应用，但仍然需要进一步的研究和改进，以提高其分解效果和精度。
DV
信号处理中，频率是信号最重要的表示。
传统的傅里叶变换分析方法并不能分析出信号的某一频率在甚么时刻出现，为此产生了能同时在时间和频率上表示信号密度和强度的时频分析，如短时傅里叶变换和小波变换等，但其基本思想都是根据傅里叶分析理论，对非线性非平稳信号的分析能力不足，受限于Heisenberg不确定原理。
HHT ( Hilbert Huang Transform)是由N. E.Huang 等人在1998 年提出的一种崭新的时频分析方法，能够对非线性非平稳的信号进行分析，同时具有良好自适应性的特点。
其本质是对信号进行平稳化处理，将具有不同时间尺度的信号逐级分解开来。
HHT 方法在各领域已得到了广泛应用，但依然存在一些不足，例如易产生虚假分量和模态混叠等。
针对传统经验模式( Empirical Mode Decomposit iON，EMD)分解方法所导致的模态混叠现象，法国以Flandrin 为首的EMD 算法研究小组和Huang 本人的研究小组通过对EMD 分解白噪声结果统计特性的大量研究，提出通过加噪声辅助分析( NADA ) 的EEMD ( EnsembleEmpirical Mode Decomposition) 方法，将白噪声加入信号来补充一些缺失的尺度，在信号分解中具有良好的表现。
EEMD仿真系统的实现利用了Matlab 平台，通过GUI 控件实现了系统设计，能直观方便地进行比较分析，验证了EEMD 在抗混叠方面较原有方法的改进。
1 经验模式分解( EMD) 和IMFHHT 方法包含两个主要步骤:( 1) 对原始数据进行经验模式分解( EMD) ，把数据分解为满足Hilbert 变换要求的n 阶本征模式函数( IMF) 和残余函数之和。
( 2) 对每一阶IMF 进行Hilbert 变换，得到瞬时频率，从而求得时频图。
函数必须关于时间轴局部对称，且其过零点与极值点个数相同。
此类函数被称为固有模态函数( Int rinsicMode Function，IMF) 。
EMD
算法概述
• EMD方法基于信号本身的局部特征时间尺度，把原始信号进行平稳化处理，将 复杂的信号分解成有限个具有不同特征尺度的数据序列，每一个序列即为一个 本征模态函数(Intrinsic Mode Function)分量，IMF反映了原始信号的本质和真实 信息。信号经EMD分解之后，其瞬时频率也具有了物理意义，因此，EMD算 法是一种非常适用于非平稳、非线性数据序列的复杂信号处理方法。
数学基础，如：正交性、收敛性、完备性、唯一性等EMD特性，试验方法求证一 些特性，而不能进行数学上的证明，甚至于至今为止都无法很好的解释“什么信 号能进行EMD分析，什么信号不能进行EMD分析”。然而对于本征模态函数， 也仅仅只能通过窄带信号的过零点与过极值点的关系以及非常有限的可用例子的 经验中获得IMF定义，其效果很难令人满意。尽管大部分的例子都表明了EMD结 果的直观合理性，但是其理论框尚待改善。
式子：
n
x (t ) c i rn
i 1
残差 r n是信号 x(t) 的集中趋势，IMFs(c1,,cn )分别包含了信号不同时间特征尺度大小的
成分，其尺度依次由小到大，因此，各分量也就相应地包含了从高到低的不同频率段的
成分。每个频率段包含的频率成分是不同的，它们随 x(t)的变化而变化。
• EMD只需要根据信号的时间特征尺度自适应的对信号进行分解。信号经EMD分解所得到的本 征模态函数均代表着信号不同尺度的特征。因为对于每个本征模态函数，连续两个极值点之间 定义了信号局部波动特征，这就反映了信号在不同尺度的特性。
本征模态函数(Intrinsic Mode Function)
一般认为，一个本征模函数IMF必须满足以下两个条件: （1）在整个信号上，极值点的个数和过零点的个数相等或至多相差一个； （2）在任意时刻，由局部极大值点和局部极小值点分别形成的上、下包络线的均值为零，也即是 说，上、下包络线相对于时间轴是局部对称的。 通常情况下，实际信号都是复杂信号，并不满足上述条件，因此，Huang进行了以下假设： （1）任何信号都是由若干本征模态函数组成的； （2）各个本征模态函数既可是线性的，也可是非线性的，各本征模态函数的局部极值点和零点相 同，同时上、下包络关于时间轴局部对称； （3）在任何时候一个信号都可以包含若干个本征模态函数，若各模态函数之间相互混叠，就组成 了复合信号。
HHT
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默认计算的信号是从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).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(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=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变换直接求该信号的瞬时频率1.clear;clc;clf;2.%假设待分析的函数是z=t^33.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(sx);17.sp=dx./dt;18.plot(t(1:N-1),sp)19.title('瞬时频率')20.复制代码小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。
5.1 EMD方法介绍及实证分析
EMD方法介绍及实证分析目录1.总体经验模式分解方法介绍 (1)1.1 EMD方法的引入 (1)1.2 EMD的基本理论和方法 (2)1.3 EEMD (3)2.实证分析 (4)2.1汇率算例分析 (4)2.2 基于EMD和GARCH模型的股价预测分析 (10)2.2.1 研究对象与数据选取 (10)2.2.2 EMD分解及分析 (10)2.2.3 自回归模型的拟合和预测 (15)2.2.4 GARCH模型的拟合和预测 (18)2.2.5 预测数据重组 (23)参考文献 (24)1.总体经验模式分解方法介绍1.1 EMD方法的引入近年来，小波变换（Wavelet Transformation , WT）理论在股票市场系统变量的多时间尺度分析与建模中取得了丰富的成果。
小波变换在时域和频域都具有良好的多分辨率分析能力，被誉为数学显微镜。
但小波变换实质上是一种窗口可调的傅立叶变换，其小波窗内的信号必须是平稳的，因而没有从根本上摆脱傅立叶分析的局限，小波变换虽然能够在频域和时域内同时得到较高的分辨率，但仍然存在一定的限制，这种限制通常会造成很多虚假的谐波，且小波基函数的选择对小波分解结果有显著的影响。
针对小波变换的不足，1998年，Huang等人提出来一种全新的多分辨率信号分析方法——经验模态分解（Empirical Mode Decomposition , EMD）。
EMD是基于信号局部特征时间尺度，从原信号中提取本征模态函数（Intrinsic Mode Function , IMF）。
在线性框架下基于EMD得到的Hilbert谱与小波谱具有相同的表现特性，而Hilbert谱在频域和时域内的分辨率都远高于小波谱，依此得到的分析结果可以更准确地反映系统原有的物理特性。
由于EMD方法比小波变换有更强的局部表现力，所以在处理非线性、非平稳信号时，EMD方法是一种更有效的方法，而金融时间序列（如股价、股价指数、收益率等）就是一类典型的非线性、非平稳时间序列。
emd经验模态分解程序
% in case the current mode is so small machine precision can cause
% spurious extrema to appear
if (max(abs(m))) < (1e-10)*(max(abs(x)))
disp(['mode ',int2str(k),', iteration ',int2str(nbit)])
if exist('s','var')
disp(['stop parameter mean value : ',num2str(s)])
end
elseif FIXE_H
stop_count = 0;
[stop_sift,moyenne] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);
else
[stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);
[indmin,indmax,indzer] = extr(y);
nem(k) = length(indmin)+length(indmax);
nzm(k) = length(indzer);
MODE_COMPLEX=0;
display_sifting=0;
mask=0;
% [x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:});
常见不同模态信号分解方法探讨
常见不同模态信号分解方法探讨邢昀;荣剑【摘要】经验模态分解(EMD)是一种自适应的信号时频分析方法,它把信号分解成一系列本征模态函数(IMF)和残差分量.集合经验模态分解方法(EEMD)是通过向原始信号中加入高斯白噪声,来抑制经验模态分解过程中存在的模态混叠现象.补充的EEMD(CEEMD)是通过向目标信号添加成对的符号相反的白噪声,来确保信号分解具有真实的物理意义.改进的集合经验模态分解(MEEMD)结合CEEMD与排列熵(PE)算法在抑制模态混叠方面取得理想的结果,并解决计算量大的问题.变分模态分解(VMD)是在EMD的基础上发展出来的一种新型信号处理方法,它进一步避免模态混叠现象并且有着更高的运算效率.讨论EMD、EEMD、CEEMD、MEEMD、VMD在信号分解处理时的效果差异.【期刊名称】《现代计算机（专业版）》【年(卷),期】2018(000)036【总页数】5页(P7-11)【关键词】经验模态分解(EMD);集合经验模态分解(EEMD);补充的集合经验模态分解(CEEMD);改进的集合经验模态分解(MEEMD);变分模态分解(VMD)【作者】邢昀;荣剑【作者单位】西南林业大学大数据与智能工程学院,昆明650224;西南林业大学大数据与智能工程学院,昆明650224【正文语种】中文0 引言在信号处理领域，从1882年傅里叶提出傅里叶级数，到1965年图基和库利发表“快速傅里叶变换算法”以来，该学科蓬勃发展。
在经典的信号处理理论中，时域和频域的关系是信号处理中的一个重要关系，傅里叶变换和傅里叶反变换在信号时域和频域之间建立起了沟通的桥梁[1]。
然而傅里叶变换只是一种全局意义上的变换，所以在分析平稳信号时候比较有效，但在实际应用中，大多数信号都是非平稳信号[2]。
非平稳信号同平稳信号相比，其分布参数或分布律随时间发生了变化。
为了处理非平稳信号，人们在傅里叶变换的基础上对其进行不断的改进和拓展，其中时频分析方法是重要分支之一。
EMD方法
经验模态分解(Empirical Mode Decomposition,简称EMD)特点：1. 适合于分析非线性、非平稳信号序列，具有很高的信噪比2. 依据数据自身的时间尺度特征来进行信号分解, 无须预先设定任何基函数，在理论上可以应用于任何类型的信号的分解EMD方法假设任何信号都由不同的本征模态函数(IMF)组成，每个IMF可以是线性的，也可以是非线性的，IMF分量必须满足下面两个条件：一是其极值点个数和过零点数相同或最多相差一个，二是其上下包络关于时间轴局部对称。
这样任何一个信号就可以分解为有限个IMF之和。
EMD分解过程基于以下假设：(1)信号最少有一个极大值和一个极小值；(2)时域特性由极值间隔决定；(3)如果数据序列完全缺乏极值但是仅包含拐点，那么它也可通过求导一次或多次来表示极值点，而最终结果可以由这些成分求积分来获得。
具体方法是由一个“筛选”过程完成的：(1)首先找出信号s(t)所有的极大值点并将其用三次样条函数拟合成原数据序列上的包络线，再找出所有的极小值点并将其用三次样条函数拟合成原数据。
(2)计算上下包络线的均值，记为m1(t)，把原数据序列s(t)减去该均值即可得到一个去掉低频的新数据序列h1:(3)因为h1(t)一般仍不是一个IMF分量序列，为此需要对它重复进行上述处理过程。
重复进行上述处理过程k次，直到h1(t)符合IMF的定义要求，所得到的均值趋于零为止，这样就得到了第1个IMF分量c1(t)，它代表信号s(t)中最高频率的分量：(4)将c1(t)从s(t)中分离出来，即得到一个去掉高频分量的差值信号,即有：将r1(t)作为原始数据，重复步骤（1）、（2）、和（3），得到第二个IMF 分量c2(t)，重复n次，得到n个IMF分量。
这样就有：当满足给定的终止条件（通常使rn(t)成为一个单调函数）时，循环结束，由上面两个式子可以得到：其中，rn(t)为残余函数，代表信号的平均趋势。
基于EMD降噪的燃烧器火焰静电信号能量熵分析
第52卷第1期2021年1月中南大学学报(自然科学版)Journal of Central South University (Science and Technology)V ol.52No.1Jan.2021基于EMD 降噪的燃烧器火焰静电信号能量熵分析李珊1，闫勇2，吴佳丽1，钱相臣1(1.华北电力大学控制与计算机工程学院，北京，102206；2.肯特大学工程与数字艺术学院，坎特伯雷，肯特CT27NT ，英国)摘要：采用灵敏度高、结构简单、环境适应性强的非侵入式静电传感器阵列测量甲烷及生物质火焰的能量变化特性。
基于经验模态分解的降噪方法对测量信号进行去噪处理，以获得信号的能量熵用于表征火焰中带电颗粒运动复杂度。
实验中对甲烷−空气预混和扩散火焰以及甲烷助燃生物质火焰进行测量。
研究结果表明：甲烷扩散火焰带电颗粒运动复杂度随燃料流量的增大而增大，相同甲烷流量下，空气流量越大预混火焰颗粒运动复杂度越高，生物质火焰中带电颗粒运动复杂度比气态甲烷火焰的高，火焰中部带电颗粒运动复杂度最高。
关键词：可再生低碳燃料；火焰检测；静电传感器阵列；经验模态分解；能量熵中图分类号：TP212文献标志码：A开放科学(资源服务)标识码(OSID)文章编号：1672-7207（2021）01-0285-09Energy entropy analysis of flame signals obtained by anelectrostatic sensor array based on EMD denoising methodLI Shan 1,YAN Yong 2,WU Jiali 1,QIAN Xiangchen 1(1.School of Control and Computer Engineering,North China Electric Power University,Beijing 102206,China;2.School of Engineering and Digital Arts,University of Kent,Kent CT27NT,U.K.)Abstract:A non-invasive electrostatic sensor array with high sensitivity,simple structure and good environmental adaptability was used to measure the energy characteristics of burner flames in different combustion tests.A denoising method based on empirical mode decomposition(EMD)was adopted to denoise the original sensor signals before calculating the EMD energy entropy of the flame,which was used to evaluate the motion complexity of charged particles in a flame.Experimental tests were conducted under various flame conditions,such as methane-air diffusion flames,methane-air premixed flames and methane-assisted biomass flames.The results show that the motion complexity of charged particles in diffusion and premixed flames increases with theDOI:10.11817/j.issn.1672-7207.2021.01.029收稿日期：2020−10−12；修回日期：2020−11−13基金项目(Foundation item)：国家自然科学基金资助项目(61673170,51827808)；中央高校基本科研业务费专项资金资助项目(2019MS023)(Projects(61673170,51827808)supported by the National Natural Science Foundation of China;Project (2019MS023)supported by the Fundamental Research Funds for the Central Universities)通信作者：钱相臣，博士，副教授，从事多相流与火焰检测技术、智能仪表与工业参数检测研究；E-mail:***************.cn引用格式：李珊,闫勇,吴佳丽,等.基于EMD 降噪的燃烧器火焰静电信号能量熵分析[J].中南大学学报(自然科学版),2021,52(1):285−293.Citation:LI Shan YAN Yong,WU Jiali.Energy entropy analysis of flame signals obtained by an electrostatic sensor array based on EMD denoising method[J].Journal of Central South University(Science and Technology),2021,52(1):285−293.第52卷中南大学学报(自然科学版)flow rate of methane.The motion complexity of charged particles in biomass flames is higher than that of methane flames.The motion complexity of charged particles in the middle region of a flame is the highest.Key words:renewable low carbon fuel;flame monitoring;electrostatic sensor arrays;empirical mode decomposition;energy entropy能源的低碳化和可再生化是全球能源结构转型的发展趋势。
EDA分解过程
3.2.1 EMD分解过程clc;t=0:0.005:10;%ScopeData.time;Delta=0.005;Fs=1/Delta; % 采样频率x=2*exp(-0.02*t).*sin(40*pi*t)+3*exp(-0.03*t).*sin(80*pi*t)+0.5*randn;%已知的系统相应函数figure;plot(t,x);xlabel('TIME/ t(s)'),ylabel('x(t)')title('原史信号图');c = x(:)'; % copy of the input signal (as a row vector)N = length(x);imf = []; % Matrix which will contain the successive IMF， and the residue while (1) % the stop criterion is tested at the end of the loop% inner loop to find each imfh = c; % at the beginning of the sifting process， h is the signalSD = 1;while SD > 0.3d = diff(h); % 微分maxmin = []; % 无区分存储极值for i=1:N-2if d(i)==0 % we are on a zeromaxmin = [maxmin, i];elseif sign(d(i))~=sign(d(i+1))% we are straddling a zero somaxmin = [maxmin, i+1]; % define zero as at i+1 (not i)endendif size(maxmin,2) < 2 % then it is the residue求数列的长度，%共有几个元素；breakend% divide maxmin into maxes and minsif h(maxmin(1))>h(maxmin(2)) % first one is a max not a minmaxes = maxmin(1:2:length(maxmin));mins = maxmin(2:2:length(maxmin));else % is the other way aroundmaxes = maxmin(2:2:length(maxmin));mins = maxmin(1:2:length(maxmin));end% make endpoints both maxes and minsm=length(maxes);n=length(mins);maxes = [1 maxes N];mins = [1 mins N];% spline interpolate to get max and min envelopes； form imf maxenv = spline(maxes,h(maxes),1:N);minenv = spline(mins, h(mins),1:N);m = (maxenv + minenv)/2; % mean of max and min enveloppesprevh = h; % copy of the previous value of h before modifying it h = h - m; % substract mean to h% calculate standard deviationeps = 0.0000001; % to avoid zero valuesSD = sum( ((prevh - h).^2) ./ (prevh.^2 + eps) );endimf = [imf; h]; % store the extracted IMF in the matrix imfif size(maxmin,2) < 2breakendc = c - h; % substract the extracted IMF from the signalend[n1,n2]=size(imf);R=[];for k=1:n1r=corrcoef(x,imf(k,:));R=[R r];endv=0;for k=1:n1u(k)=R(1,2*k);if abs(u(k))>=abs(u(1))v=abs(u(k));endendw=10;p=v/w;E=0;X=[];S=[];for k=1:n1-1if abs(R(1,2*k))>pE=E+1;X=[X;imf(k,:)];figure;plot(t,imf(k,:)) %IMF分量的时域图xlabel('TIMF/ t(s)');ylabel('imf');np=1001; %输出数据长度t1=0:1/Fs:(np-1)/Fs; %建立离散输出时间向量x1=imf(k,:);nt=length(x1); %取输入数据长度s=1.5*std(x1); %设置截取振幅为输入信号标准差的1.5倍 %获取输入信号的子样本函数进行叠加m=0;y=zeros(1,np);for i=2:nt-npa=abs(x1(i-1)-s);b=abs(x1(i)-s);c=abs(x1(i+1)-s);if b<a & b<cy(1:np)=y(1:np)+x1(i:i+np-1);m=m+1;endendendendsize(S),for k=1:n1-1if abs(R(1,2*k))>=pfigure;nfft=512;Y=fft(imf(k,:),nfft);f=Fs*(0:(nfft/2-1))/nfft;fftY=abs(Y)/nfft;plot(f,fftY(1:nfft/2)),grid %IMF分量的频域图 xlabel('Frequence/Hz');ylabel('Amp');elseimf(n1,:)=imf(n1,:)+imf(k,:);endendimf=X;save imf;3.2.1 Hilbert变换求模态参数load imf; %load前面EMD分解得到的IMF数据Delta=0.005;Fs=1/Delta; % 采样频率a=imf(1,:);plot(a);ya=hilbert(a');p=gradient(a,0.005);p1=gradient(p,0.005);Ya=hilbert(a);A=abs(Ya);A1=gradient(A,0.005);Yx=hilbert(x);yx=imag(Yx);ya1=gradient(ya,0.005);ya2=gradient(ya1,0.005);wr=instfreq(ya);figure;Wr=wr*1256; % 瞬时频率求解,instfreq函数算出的频率应乘以2*pi和%采样频率才可得出圆频率np=length(Wr);t1=0:1/Fs:(np-1)/Fs;plot(t1(100:900),Wr(100:900)/2/pi);hold onxlabel('时间(s)');ylabel('频率（Hz）');title('频率实测图');axis([0.5,5,0,50]); hold on;a=imf(2,:);ya=hilbert(a');p=gradient(a,0.005);p1=gradient(p,0.005);Ya=hilbert(a);A=abs(Ya);A1=gradient(A,0.005);Yx=hilbert(x);yx=imag(Yx);ya1=gradient(ya,0.005);ya2=gradient(ya1,0.005);wr=instfreq(ya);Wr=wr*1256; % 瞬时频率求解,instfreq函数算出的频率应乘以2*pi和%采样频率才可得出圆频率np=length(Wr);t1=0:1/Fs:(np-1)/Fs;plot(t1(100:900),Wr(100:900)/2/pi);hold onxlabel('时间(s)');ylabel('频率（Hz）');title('频率实测图');function [fnormhat,t]=instfreq(x,t,L,trace);%INSTFREQ Instantaneous frequency estimation.% [FNORMHAT,T]=INSTFREQ(X,T,L,TRACE) computes the instantaneous % frequency of the analytic signal X at time instant(s) T, using the% trapezoidal integration rule.% The result FNORMHAT lies between 0.0 and 0.5.%% X : Analytic signal to be analyzed.% T : Time instants (default : 2:length(X)-1).% L : If L=1, computes the (normalized) instantaneous frequency% of the signal X defined as angle(X(T+1)*conj(X(T-1)) ;% if L>1, computes a Maximum Likelihood estimation of the% instantaneous frequency of the deterministic part of the signal% blurried in a white gaussian noise.% L must be an integer (default : 1).% TRACE : if nonzero, the progression of the algorithm is shown% (default : 0).% FNORMHAT : Output (normalized) instantaneous frequency.% T : Time instants.%% Examples :% x=fmsin(70,0.05,0.35,25); [instf,t]=instfreq(x); plot(t,instf)% N=64; SNR=10.0; L=4; t=L+1:N-L; x=fmsin(N,0.05,0.35,40);% sig=sigmerge(x,hilbert(randn(N,1)),SNR);% plotifl(t,[instfreq(sig,t,L),instfreq(x,t)]); grid;% title ('theoretical and estimated instantaneous frequencies');%% See also KAYTTH, SGRPDLAY.% F. Auger, March 1994, July 1995.% Copyright (c) 1996 by CNRS (France).%% ------------------- CONFIDENTIAL PROGRAM --------------------% This program can not be used without the authorization of its% author(s). For any comment or bug report, please send e-mail to% f.auger@if (nargin == 0),error('At least one parameter required');end;[xrow,xcol] = size(x);if (xcol~=1),error('X must have only one column');endif (nargin == 1),t=2:xrow-1; L=1; trace=0.0;elseif (nargin == 2),L = 1; trace=0.0;elseif (nargin == 3),trace=0.0;end;if L<1,error('L must be >=1');end[trow,tcol] = size(t);if (trow~=1),error('T must have only one row');end;if (L==1),if any(t==1)|any(t==xrow),error('T can not be equal to 1 neither to the last element of X'); elsefnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi); end;elseH=kaytth(L);if any(t<=L)|any(t+L>xrow),error('The relation L<T<=length(X)-L must be satisfied'); elsefor icol=1:tcol,if trace, disprog(icol,tcol,10); end;ti = t(icol); tau = 0:L;R = x(ti+tau).*conj(x(ti-tau));M4 = R(2:L+1).*conj(R(1:L));diff=2e-6;tetapred = H * (unwrap(angle(-M4))+pi);while tetapred<0.0 , tetapred=tetapred+(2*pi); end;while tetapred>2*pi, tetapred=tetapred-(2*pi); end;iter = 1;while (diff > 1e-6)&(iter<50),M4bis=M4 .* exp(-j*2.0*tetapred);teta = H * (unwrap(angle(M4bis))+2.0*tetapred);while teta<0.0 , teta=(2*pi)+teta; end;while teta>2*pi, teta=teta-(2*pi); end;diff=abs(teta-tetapred);tetapred=teta; iter=iter+1;end;fnormhat(icol,1)=teta/(2*pi);end;end;end;。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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默认计算的信号是从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).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(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=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变换直接求该信号的瞬时频率1.clear;clc;clf;2.%假设待分析的函数是z=t^33.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(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*10*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,:).^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',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,'fontsize',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);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.ylabel(['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);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.figure(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.xlabel('频率 / 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.151.if min(size(x)) == 1152.if size(x,2) == 1153.x = x';154.if nargin < 2155.t = 1:size(x,2);156.end157.end158.Nmodes = 1;159.else160.Nmodes = size(x,1);161.end162.163.lt=length(t);164.165.tt=t((l+1):(lt-l));166.167.for i=1:Nmodes168.169.an(i,:)=hilbert(x(i,:)')';170.f(i,:)=instfreq(an(i,:)',tt,l)'; 171.A=abs(an(:,l+1:end-l));172.173.if aff174.disprog(i,Nmodes,max(Nmodes,100))175.end176.177.end178.179.180.function disp_hhs(im,t,inf)181.182.% DISP_HHS(im,t,inf)183.% displays in a new figure the spectrum contained in matrix "im" 184.% (amplitudes in log).185.%186.% inputs : - im : image matrix (e.g., output of "toimage") 187.% - t (optional) : time instants (e.g., output of "toimage") 188.% - inf (optional) : -dynamic range in dB (wrt max)189.% default : inf = -20190.%191.% utilisation : disp_hhs(im) ; disp_hhs(im,t) ; disp_hhs(im,inf) 192.% disp_hhs(im,t,inf)193.194.figure195.colormap(bone)196.colormap(1-colormap);197.198.if nargin==1199.inf=-20;200.t = 1:size(im,2);201.202.end203.204.if nargin == 2205.if length(t) == 1206.inf = t;207.t = 1:size(im,2);208.else209.inf = -20;210.end211.end212.213.if inf >= 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.ylabel(['normalized frequency'])228.title('Hilbert-Huang spectrum')229.function [f,z]=fftfenxi(t,y)230.L=length(t);N=2^nextpow2(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*f)得到频移后[-2,2]之间的频谱值237.Y=fft(y);238.z=sqrt(Y.*conj(Y));239.240.241.242.复制代码4.总结。