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算法的可靠性和稳定性。
总之，经验模态分解算法是一种有效的信号分解方法，能够提供信号的局部特征表示。
它在很多领域有广泛的应用，但仍然需要进一步的研究和改进，以提高其分解效果和精度。
- 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.总结。