希尔伯特-黄变换说明及程序(标准程序)

Hilbert-Huang变换1.来自wikipedia的解释希尔伯特-黄转换维基百科，自由的百科全书跳转到：导航,搜索希尔伯特-黄转换（Hilbert-HuangTransform），由台湾中央研究院院士黄锷（NordenE.Huang）等人提出，将欲分析资料分解为本质模态函数（intrinsicmodefunctions,IMF），这样的分解流程称为经验模态分解(EmpiricalModeDecomposition,EMD)的方法。然后将IMF作希尔伯特转换(HilbertTransform)，正确地获得资料的瞬时频率。此方法处理对象乃针对非稳态与非线性讯号。与其他数学转换运算（如傅立叶变换）不同，希尔伯特-黄转换算是一种应用在数据资料上的算法，而非理论工具。目录1本质模态函数(IMF)2经验模态分解(EMD)3结论4相关条目5参考文献6外部链接[编辑]本质模态函数(IMF)任何一个资料，满足下列两个条件即可称作本质模态函数。⒈局部极大值(localmaxima)以及局部极小值(localminima)的数目之和必须与零交越点(zerocrossing)的数目相等或是最多只能差1，也就是说一个极值后面必需马上接一个零交越点。⒉在任何时间点，局部最大值所定义的上包络线(upperenvelope)与局部极小值所定义的下包络线，取平均要接近为零。因此，一个函数若属于IMF，代表其波形局部对称于零平均值。此类函数类似于弦波（sinusoid-like），但是这些类似于弦波的部分其周期与振幅可以不是固定。因为，可以直接使用希尔伯特转换，求得有意义的瞬时频率。[编辑]经验模态分解(EMD)EMD算法流程图建立IMF是为了满足希尔伯特转换对于瞬时频率的限制条件之前置处理，也是一种转换的过程。我们将IMF来做希尔伯特转换可以得到较良好的特性，不幸的是大部分的资料并不是IMF，而是由许多弦波所合成的一个组合。如此一来，希尔伯特转换并不能得到正确的瞬时频率，我们便无法准确的分析资料。为了解决非线性（non-linear）与非稳态（non-stationary）资料在分解成IMF时所遇到的困难，便发展出EMD。经验模态分解是将讯号分解成IMF的组合。经验模态分解是借着不断重复的筛选程序来逐步找出IMF。以讯号为例，筛选程序的流程概述如下:步骤1:找出中的所有局部极大值以及局部极小值，接着利用三次样条(cubicspline)，分别将局部极大值串连成上包络线与局部极小值串连成下包络线。步骤2:求出上下包络线之平均，得到均值包络线。步骤3:原始信号与均值包络线相减，得到第一个分量。步骤4:检查是否符合IMF的条件。如果不符合，则回到步骤1并且将当作原始讯号，进行第二次的筛选。亦即重复筛选次直到符合IMF的条件，即得到第一个IMF分量，亦即步骤5:原始讯号减去可得到剩余量，表示如下式步骤6:将当作新的资料，重新执行步骤1至步骤5，得到新的剩余量。如此重复次...当第个剩余量已成为单调函数(monotonicfunction)，将无法再分解IMF时，整个EMD的分解过程完成。原始讯号可以表示成个IMF分量与一个平均趋势(meantrend)分量的组合，亦即如此一来，原始资料便分解成n个IMF和一个趋势函数，我们便可将IMF做希尔伯特转换来进行瞬时频率的分析。[编辑]结论傅立叶变换是将一个讯号分解成无限多个弦波来分析资料，但是希尔伯特-黄转换则是将一个讯号分解成数个近似于弦波的讯号（周期、振幅不固定）和一个趋势函数来做分析。两者各有其优缺点，整理如下优点：1.避免复杂的数学运算2.可分析频率会随时间变化的讯号3.较适于分析气候、经济等具有趋势的资料4.可以找出一个函数的趋势缺点：1.缺乏严谨的物理意义2.需要复杂的递回，运算时间反而比短时距傅立叶变换要长3.希尔伯特转换未必能正确计算出本质模态函数之瞬时频率4.无法使用快速傅立叶变换5.只有在特例（组合较简单的资料）时使用希尔伯特-黄转换较快[编辑]相关条目希尔伯特转换包络线样条时频分析[编辑]参考文献NordenE.Huang,etal.TheempiricalmodedecompositionandtheHilbertspectrumfornonlinearandnon-stationarytimeseriesanalysis.Proc.R.Soc.Lond.A(1998)454,903–995(Link)陈韦佑,以希尔伯特-黄转换法为GPS接收机抑制调频干扰,国立台湾大学电机工程研究所硕士论文,2007.Jian-JiunDing,Timefrequencyanalysisandwavelettransformclassppt,theDepartmentofElectricalEngineering,NationalTaiwanUniversity(NTU),Taipei,Taiwan,2009.Jian-JiunDing,时频分析与小波转换,availablein[编辑]外部链接（正体中文）国立中央大学数据分析方法研究中心（正体中文）HHT标准程序下载2.Matlab代码分解代码（emd.m）functionimf=emd(x)%EmpiricialModeDecomposition(Hilbert-HuangTransform)%imf=emd(x)%Func:findpeaksx=transpose(x(:));imf=[];while~ismonotonic(x)x1=x;sd=Inf;while(sd0.1)|~isimf(x1)s1=getspline(x1);s2=-getspline(-x1);x2=x1-(s1+s2)/2;sd=sum((x1-x2).^2)/sum(x1.^2);x1=x2;endimf{end+1}=x1;x=x-x1;endimf{end+1}=x;%FUNCTIONSfunctionu=ismonotonic(x)u1=length(findpeaks(x))*length(findpeaks(-x));ifu10,u=0;else,u=1;endfunctionu=isimf(x)N=length(x);u1=sum(x(1:N-1).*x(2:N)0);u2=length(findpeaks(x))+length(findpeaks(-x));ifabs(u1-u2)1,u=0;else,u=1;endfunctions=getspline(x)N=length(x);p=findpeaks(x);s=spline([0pN+1],[0x(p)0],1:N);2.2找极值代码（findpeaks.m）functionn=findpeaks(x)%Findpeaks.%n=findpeaks(x)n=find(diff(diff(x)0)0);u=find(x(n+1)x(n));n(u)=n(u)+1;2.3绘制时-频曲线以及尺度分解代码（plot_hht.m）functionplot_hht(x,Ts)%PlottheHHT.%plot_hht(x,Ts)%Ts：timeinterval（sec）%::Syntax%ThearrayxistheinputsignalandTsisthesamplingperiod.%Func:emd%GetHHT.imf=emd(x);fork=1:length(imf)b(k)=sum(imf{k}.*imf{k});th=angle(hilbert(imf{k}));d{k}=diff(th)/Ts/(2*pi);end[u,v]=sort(-b);b=1-b/max(b);%Settime-frequencyplots.N=length(x);c=linspace(0,(N-2)*Ts,N-1);fork=v(1:2)figure,plot(c,d{k},'k.','Color',b([kkk]),'MarkerSize',3);set(gca,'FontSize',8,'XLim',[0c(end)],'YLim',[01/2/Ts]);xlabel('Time'),ylabel('Frequency');end%SetIMFplots.M=length(imf);N=length(x);c=linspace(0,(N-1)*Ts,N);fork1=0:4:M-1figurefork2=1:min(4,M-k1),subplot(4,1,k2),plot(c,imf{k1+k2});set(gca,'FontSize',8,'XLim',[0c(end)]);endxlabel('Time');end