HHT变换讲义EMD讲解

1、HHT变换讲义1.1简介传统的信号处理方法，如傅立叶分析是一种纯频域的分析方法。它用频率不 同的各复正弦分量的叠加来拟合原函数，也即用F 来分辨f 。而F 在有限频域上的信息不足以确定在任意小范围内的函数f ，特别是非平稳信号在时间轴上的任何突变，其频谱将散布在整个频率轴上。 而且，非平稳动态信号的 统计特性与时间有关，对非平稳信号的处理需要进行时频分析，希望得到时域和 频域中非平稳信号的全貌和局域化结果。在傅立叶变换中，人们若想得到信号的 时域信息，就得不到频域信息。反之亦然。后来出现的小波（Wavelet）变换通过一种可伸缩和平移的小波对信号变换，从而达到时频局域化分析的目的。但这种变换实

2、际上没有完全摆脱傅立叶变换的局限，它是一种窗口可调的傅立叶变 换，其窗内的信号必须是平稳的。另外，小波变换是非适应性的，小波基一旦选 定，在整个信号分析过程中就只能使用这一个小波基了。HHT （ Hilbert-Hua ng Tran sform）技术是（1998年由 NASA 的 Norden E Hua ng 等提出的新的信号处理方法。该方法适用于非线性非平稳的信号分析， 被认为是 近年来对以傅立叶变换为基础的线性和稳态谱分析的一个重大突破。目前 HHT 技术已用于地球物理学和生物医学等领域的研究，并取得了较好的结果。存在的问题尽管HHT技术在处理非线性、非稳态信号方面有很大的优势，但是这

3、个方法 本身还是有许多的问题有待进一步研究。正如Hua ng在文章中指出的那样，对于这种新的信号处理方法，其基的完 备性还需要严密的证明。另外，在做Hilbert变换时出现的边界效应也需要更好的 方法来解决。但是，HHT技术中最严重，也是现今研究的最多的是 EMD分解中 的包络过程。从对EMD分解方法的介绍可以看出，包络线的构造影响着整个分 解的结果，也决定了后面的Hilbert变换。Hua ng采用的三次样条插值来拟和包络线。在实际应用中，发现这样做会产生严重的边界效应，污染了原始数据。特别 是对短数据而言，这种影响使分析所得的结果失去了原有的意义。对此，Hua ng等采用的是根据信号端点处

4、的振幅和频率，分别增加两组特征波的方式进行数据延拓。Hua ng的这种延拓方法已经向NASA申请了专利。除此之外，还有人提出 了其它方法进行端点延拓。比如国家海洋局的黄大吉等提出的镜像闭合延拓法和 极值延拓法。其中镜像闭合延拓法是根据信号的分布特性，把镜子放在具有对称性的极值位置，通过镜像法把镜内信号映射成一个周期性的环形信号，不存在端 点，从根本上避免了端点效应。而极值延拓法具有镜像闭合延拓法相当的效果， 而不增加信号序列本身的长度，计算较快，尤其在处理非对称波形信号时比较优 越。还有青岛海洋大学的邓拥军等，利用神经网络分析方法来延拓数据端。这些方法对边界效应的抑制都有一定的作用，但是也都需

5、要更深一步的研究。1.2瞬时频率频率是个极其重要的物理量，定义为信号周期倒数，其物理含义显而易见。对于正弦信号，它的频率为恒值。但是对于大部分信号，它的频率是随时间变化 的函数，故提出瞬时频率概念。瞬时频率即表征信号在局部时间点上瞬态频率特 性，整个持续期上的瞬时频率反映了信号频率的时变规律。对于随机时间序列X(t)，对其进行Hilbert变换跑，可以得到丫(t)如下：Y t 二丄 PV( . )(1-1)t -1 其中，PV为柯西主值(Cauchy principal value)。该式表示Y(t)是X(t)与一的卷积。通过这个定义，X(t)和Y(t)组成了一个共轭复数对，于是可以得到一个解

6、 析信号Z(t)如下Z(t )=X(t )+iY(t 戶 a(t)eW)(1-2)a t；LX2 t Y2 t 12其中/(1-3)日(t)=arctan iiX(t)丿从理论上讲，虚部的定义方法有很多种。但是Hilbert变换为其提供了一个 唯一的虚部值，这就使得其结果成为一个解析函数。 得到了相位，就可以得到瞬 时频率，因为瞬时频率就是相位导数。=-(1-4)dt从本质上说，公式(1-1) 将 Hilbert变换定义为X(t)与丄的卷积，因此它强 调了 X(t)的局部特性。公式(1-2)的极坐标表达式更进一步澄清了这个表达式的 局部特性：它的幅值和相位的最佳局部匹配随三角函数变化。1.3

7、固有模态函数 IMF (Intrinsic Mode Function)由瞬时频率的物理意义可知，并不是任意的信号都能用瞬时频率来讨论。 只 有当信号满足只包括一种振动模式，而没有复杂叠加波的情况时才行 。实际上， 定义一个有意义的瞬时频率的必要条件就是要求函数关于局部零平均值对称，并 且零交叉点和极值点数量相同。基于此种原因，提出了固有模态函数的概念。固 有模态函数满足以下两个条件：(1) 整个数据范围内，极值点和过零点的数量相等或者相差一个；(2) 在任意点处，所有极大值点形成的包络线和所有极小值点形成的包络线 的平均值为零。第一个条件是显而易见的，它类似于平稳过程中传统的稳定且满足高斯分

8、布 的窄带信号条件。第二个条件把传统的全局条件调整到局部情况。只有满足了这个条件，得到的瞬时频率才不会因为不对称波形的存在而引起不规则波动。所以这一点是得到正确瞬时频率的必要条件。 而这一点是必须的，因为这样瞬时频率 就可以不包含由于不对称波形造成的波动。为了使用瞬时频率定义，必须要把随机数据归结为IMF组件，这样才可以为 每个IMF组件定义瞬时频率。为了将数据归结为所需的IMF组件，接下来引入经 验模式分解方法。1.4 经验模式分解 EMD ( Empirical Mode Decomposition经验模态分解方法的大体思路是利用时间序列上下包络的平均值确定“瞬时平衡位置”，进而提取固有模

9、态函数39。这种方法基于如下假设：(1) 信号至少有两个极点一-个极大值和一个极小值；(2) 信号特征时间尺度是由极值间的时间间隔来确定的；(3) 如果数据没有极值而仅有拐点，可以通过微分、分解、再积分的方法 获得IMF。在此假设基础上，Hua ng等人进一步指出：可以用经验模态分解方法将信号 的固有模态筛选出来。经验模式分解过程就是个筛选过程，实现振动模式的提取 该方法的基本思路是用波动上、下包络的平均值去确定“瞬时平衡位置”，进而 提取出固有模态函数。上、下包络线是由三次样条函数对极大值点和极小值点进 行拟合得到的。经验模式分解过程的基本过程可概括如下：(1) 寻找信号x(t)所有局部极大

10、值和局部极小值，为更好保留原序列的特性，局部极大值定义为时间序列中的某个时刻的值，它只要满足既大于前一时刻的值也大于后一时刻的值即可。局部极小值的提取同理，即该时刻的值满足既小于前 一时刻的值也小于后一时刻的值。使用三次样条函数进行拟合，获得上包络线 Xmax(t)和下包络线Xmin(t)；(2) 计算上、下包络线的均值m(t)=xmax(t)+xmin(t)/2；(3) 用原信号x(t)减去均值m(t)，得到第一个组件h(t)=x(t)-m(t)；由于原始序列的差异，组件h(t)不一定就是一个IMF，如果h(t)不满足固有模态函数 两个条件，就把h(t)当成原始信号，重复(1)-(3)，直到

11、满足条件为止，这时满足 固有模态函数条件的h(t)作为一个IMF，令Ii(t)=h(t)，至此第一个IMF已经成功 的提取了。由于剩余的r(t)=x(t)-I i(t)仍然包含具有更长周期组件的信息，因此可 以把它看成新的信号，重复上述过程，依次得到第二个I2(t)，第三个b(t),，当r(t)满足单调序列或常值序列条件时，终止筛选过程，可以认为完成了提取固 有模态函数的任务，最后的r(t)称为余项，它是原始信号的趋势项。由此可得x(t)n的表达式X t =7 Ij t r t，即原始序列是由n个IMF与一个趋势项组成。如 i壬上所述，整个过程就像筛选过程，根据时间特性把固有模态函数从信号中提

12、取出 来。1.5端点延拓HHT方法的分析质量很大程度上取决于 EMD解的质量，而在应用EMM法 时的一个非常棘手的问题是，由于信号两端不可能同时处于极大值和极小值，这样，“筛”过程中构成上下包络的三次样条函数在数据序列的两端就会出现发散 现象(如图1.1(a)所示，原始信号是中心频率为 1KHz 5个波峰、以100个数 据描绘一周期的窄带信号)，并且这种发散的结果会随着“筛”过程的不断进行 逐渐向内“污染”整个数据序列而使所得结果严重失真(如图 1.1(b)所示)。另外，在进行Hilbert变换时，信号的两端也会出现严重的端点效应。对于一个较长的数据序列来讲，可以根据极值点的情况不断抛弃两端的

13、数据来保证所 得到的包络的失真度达到最小。但对于一个数据点数少的序列来讲，这样的操作就变得完全不可行。因此，必须对信号或其极值向外进行延拓，以确保包络线抵达端点端点延拓的目的是确保上、下包络都与端点相交，以便有与每一个信号点 相对应的局部平均值。而上、下包络是由极大值和极小值连结而成的，因此只要对极大值和极小值进行延拓，而不必对信号本身进行延拓。极大值和极小值是相 间分布的，同时考虑到样条插值的要求，所以只要在信号左、右两端分别延拓两 个极大值和两个极小值即可。由于端点以外没有信号，任何延拓都是人为的，通 过对不同形式的信号的多次试验，发现一种端点延拓的方法以端点的一个特征波 为依据进行延拓具

14、有较好的效果，具体做法如下40:-ia103030 AD 5D 6070 BQ 90100乘样点萌图3.1未作端点延拓的分析效果1-5 IIIIIII1口1QS 口一启科10(a)拟合结果D 占 11.523.533.544.5051152 空牟 335 a d(b) EMD分解QD5 11.522.533.540=5115 J 3$335 A 2Tirncifrna)3-SO(a)拟合结果(b) EMD分解图3.2以端点的一个特征波为依据进行延拓的分析效果设离散信号：t t(1),t(2),t(n)电,t2，,tn(1-5)X(t) X(tJ,X(t2),X(tn) “NX, ,Xn(1-6

15、)其采样步长为-=t , X(t)有M个极大值和N个极小值，对应的序列下标(lm, ln)、 时间(Tm , Tn)和函数值(U ， V )记为：Im 珂U1),U2),IH,Im(M)(1-7)In 二I,1),In (2)Mn(N)(1-8)Tmd，U(i)=Xi =1,|,M(1-9)Tn(i)=tin， V(i) =Xin, i =1, ,N(1-10)具体作法如下：1. 左端信号左端第一个特征波包含的信号点数为k1Im(2) -， 当 lm(1)n(1) ln(2) -ln(1，当 Im In(1-11)2|lm(1)n(1)|，当 M=N=1向外延拓的两个极值的位置（Tm, Tn

16、）和数值（U , V ）为：0（0）二七：t,U(0)=U(1)(1-12)0(-1)=0(1)-2心龙U(-1)=U(1)(1-13)(0) =Tn(1)-k1：tV(0) =V(1)(1-14)人（-1）=（1）-23V(-1) = V(1)(1-15)2. 右端信号右端第一个特征波包含的信号点数为k2(1-16)h(M)lm(M -1),当 lm(M)ln(N)k2 才 ln(N)-ln(N -1),当 lm(M)Cln(N)2|lm(M)-ln(N)|， 当 M=N=1向外延拓的两个极值的位置（Tm, Tn ）和数值（U , V ）为:Tm(M 1)=Tm(M) k2 ：t,U (M

17、1) =U(M)(1-17)Tm(M 2)=Tm(M) 2k2：t,U(M 2) =U(M)(1-18)Tn(N 1) =Tn(N) k2 ：t,V(N 1)=V(N)(1-19)Tn(N 2) =Tn(N) 2k2：t,V(N 2) =V(N)(1-20)3.当端点的数值比近端点的第一个极大值大或极小值小时，要进行特殊的处理，以避免信号落到包络线之外，Tm(O)丸，U(0) 2,当 X1U(1)(1-21)Tn(O) 7,V(0)=冷当 X :V(1)(1-22)Tm(M 1)=tn, U(M 1) = X“当 Xn U(M)(1-23)Tn(N 1)=tn, V(N 1) = Xn 当 X

18、n V(N)(1-24)由图1.2可见，以端点的一个特征波为依据进行延拓， 分别在信号两端增加 两个极大值和两个极小值，从而使原始信号被延长的包络线所限制，有效地抑制 了端点效应，使得EMD解得到合理的各个IMF模态。另外，国家海洋局的黄大吉等还提出了镜像闭合延拓法，即根据信号的分布特性，把镜子放在具有对称性的极值位置， 通过镜像法把境内信号映射成一个周 期性的环形信号。这样做就不存在端点，从根本上避免了端点效应40。青岛海洋 大学的邓拥军等人利用神经网络分析的方法来进行端点延拓，在数据的两端各得 到两个附加的极大值点和两个附加的极小值点，由此有效的抑制了端点效应41。1.6 EMD结束准则根

19、据理论分析IMF必须满足两个条件，但是在实际筛选过程中，严格满足这 两个条件的信号几乎不存在，因此如果以这两个条件为依据判定IMF，可能根本就得不到结果或者要以冗长的程序执行时间为代价。因此筛选过程要小心进行，为了保证IMF组件的调幅和调频都具有物理意义,习惯上，可以通过标准同时考虑到程序可行性，必须制定一个结束筛选的准则 偏差SD来完成，它可以由连续两个筛选结果得到18:(1-25)hi(k-i)(t)和hik(t)分别表示两个连续的筛选结果，一般来说，SD值越小，所得的IMF分量的线性和稳定性就越好，典型的SD值为0.2-0.3，这样对于不同筛选过程就 有严格的SD限制了。但是在筛选过程中

20、发现，SD的这个范围过于死板，对很多 信号而言，以0.2-0.3为限制的SD值并不科学，因此在实际筛选过程中，可以以 0.2-0.3为参考值，根据实际情况进行适当调整。另外，如果不是进行信号精确 处理，可以考虑控制筛选次数来决定筛选结果，这样不仅保障了程序的有效性， 也控制了程序的执行时间，在具体问题中取得了很好的效果。1.7 Hilbert 谱完成了经验模式分解过程就得到了所有可提取的固有模态函数，在Hilbert变换的基础上，只要根据公式(1-1)-(1-4)计算瞬时频率就可以了。对固有模态 函数进行Hilbert变换后，可以用下面方式表示信号nX t 二 at exp i 打 t dt(

21、1-26)jm在这儿没有考虑余项rn，因为它只是单调函数或常量。尽管Hilbert变换可以把 单调函数看成是振动的一部分，但是其余项中的能量很小，一般情况可以不用考由希尔伯特变换得出的振幅和频率都是时间的函数， 如果用三维图形表达幅 值、频率和时间之间的关系，或者把振幅用灰度的形式显示在频率一时间平面上， 就可以得到Hilbert谱H (w,t)。如果把H (w,t)对时间积分，就得到希尔伯特边际谱 h(w):Th(w) H(w,t)dt(1-27)边际谱提供了对每个频率的总振幅的量测，表达了整个时间长度内累积的振幅。另外，作为希尔伯特边际谱的附加结果，可以得到如式(1-28 )定义的希尔 伯

22、特瞬时能量：2IE(t) H (w,t)dw(1-28)w瞬时能量提供了信号能量随时间的变换情况。事实上，如果振幅的平方对时 间积分，可以得到希尔伯特能量谱：T 2ES(w)二 o H2(w,t)dt(1-29)希尔伯特能量谱提供了对于每个频率的能量的量测，表达了每个频率在整个时间长度内所累积的能量。1.8 HHT方法的软件实现HH变换的整个过程包括EM与Hilbert变换，整个过程是通过matlab软件编 写的M文件实现的，其程序流程图如图1.3所示。、开始1.9 HHT方法的优越性1、HH方法是一种全新的信号分析方法，能够描绘出信号的时频谱图、边际 谱，能量谱等，是一种更具有适应性的时频局

23、域化分析方法28 0而传统的信号处 理方法，如傅立叶变换是一种纯时域的分析方法，它用频率从零到无穷大的各复正弦分量的叠加来拟合原函数f(t)在每个时刻的值，也即用F()在有限频域上 的信息不足以确定在很小范围内的函数f(t)，特别是非平稳信号在时间轴上的任 何突变，其频谱将散布在整个频率轴上39 0所以，这种分析方法适用于确定性的平稳信号，而在非线性、非平稳过程的处理上，傅立叶变换将只作为一种数学变化手段。而且，非平稳信号的统计特性与时间有关，对非平稳信号的处理需要进 行时频分析，希望得到时域和频域中非平稳信号的全貌和局域化结果。但在傅立 叶变换中，若想得到信号的时域信息，就得不到频域信息，反之亦然。后来出现 的小波变换通过一种可伸缩和平移的基小波对信号变换，从而达到时域局域化分析的目的，但这种变化实际上并没有完全摆脱傅立叶变换的局限，它是一种窗口可调的傅立叶变换，其窗内的信号必须是平稳的12。2、HH局部性能良好而且是自适应的，对平稳信号和非平稳信号都能进行分析。它没有固定的先验基底，分解完全基于数据本身进行。固有模态函数是基于 序列数据的时间特征尺度得出的不同的数据有不同的固有模态函数，每个固有模态函数可以认为是信号中固有的一个模态，所以通过Hilbert变换得到的瞬时频率具