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）在任何时候一个信号都可以包含若干个本征模态函数，若各模态函数之间相互混叠，就组成 了复合信号。
- 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.总结。