基于EMD的信号瞬时频率估计_刘小丹

第32卷第1期2009年3月 辽宁师范大学学报(自然科学版)Journal of Liaoning Normal University (Natural Science Edition ) Vol.32　No.1Mar.　2009 文章编号:100021735(2009)0120051207基于EMD 的信号瞬时频率估计刘小丹,　孙晓奇,　沈　滨(辽宁师范大学计算机与信息技术学院,辽宁大连　116029)收稿日期:2008209224基金项目:辽宁省教育厅科学技术研究项目(20060466)作者简介:刘小丹(19572),男,吉林蛟河人,辽宁师范大学教授,硕士.E 2mail :xdliu @摘　要:分析了信号瞬时频率的定义及其两种主要的获得信号相位的方法:解析信号法和正交模型法.提出了一种基于经验模式分解的新的瞬时频率估计方法———正交包络法.该方法计算简单,克服了正交模型法无法由一个时间函数确定两个时间函数的困难.与Hilbert 变换方法相比,正交包络法使边界问题得到了明显改善.实验证明这是一种有效的瞬时频率估计方法.关键词:瞬时频率;正交包络法;EMD ;Hilbert 变换中图分类号:TP202.4 文献标识码:A根据Fo urier 分析理论,任何一个平稳信号都可以表示为多个谐波的加权和,对于谐波的某一特定频率,其幅值和相位是常数.而对于非平稳信号,由于其谱特性是随时间变化的,因此不能简单地用Fourier 变换作为非平稳信号的分析工具[1],平稳信号的频率概念也就无法准确解释非平稳信号的时变特性,于是就需要引入一个随时间变化的频率的概念,即瞬时频率.瞬时频率的一个重要特性是作为时间的函数,用它可以确定信号谱峰的位置.基于这一特性,瞬时频率的概念有着极其重要的应用,因此瞬时频率的估计也就成为许多实际的信号处理应用中一项很有意义的工作.一些信息探测系统只要系统与目标之间有相对运动,多普勒效应就会使频率改变,传播媒质的扰动也会使频率变化,雷达、声呐、移动通信、医疗设备和天文观测都存在这一问题.以雷达信号处理为例,其主要目的是对目标实行检测、跟踪和成像,而像军用飞机一类的目标为了逃避被跟踪,其径向速度是随时间改变的,这使得雷达的多普勒频率具有非平稳的谱.因此,跟踪这类目标需要用到瞬时频率估计技术.瞬时频率估计技术也应用于生物医学.例如,血流的多普勒变化直接关系到心脑血管疾病的诊断.同时,在地震信号处理中,可以利用瞬时频率来确定不同的地质构造.在语音处理等其他诸多领域都有瞬时频率估计技术的应用,详见文献[223].从物理学的角度,信号可以分为单分量信号和多分量信号.单分量信号在任意时刻都只有一个频率,该频率称为信号的瞬时频率,而多分量信号则在某些时刻具有多个不同的瞬时频率.瞬时频率的定义最早是由Carson 和Fry 在研究调频信号时分别提出的,在Gabor 提出了解析信号的概念之后,Ville 将二者结合起来,提出了现在普遍接受的实信号的瞬时频率的定义[4],即:实信号的瞬时频率就是该信号所对应的解析信号的相位关于时间的导数.上述定义只对单分量信号有意义.下面分析一下将瞬时频率定义为复信号相位关于时间的导数的原因.设一复信号c (t )=A (t )e j φ(t ),A (t )、φ(t )分别称为信号c (t )的幅度和相位.c (t )的频谱为C (ω)=12π∫+∞-∞c (t )e -j ωt d t c (t )的总能量E =∫+∞-∞|c (t )|2d t =∫+∞-∞|C (ω)|2d ω 于是,归一化的函数|c (t )|2/E 和|C (ω)|2/E 可分别作为信号c (t )在时域和频域的能量密度函数,从而得到信号频谱C (ω)的平均频率: 〈ω〉=1E ∫+∞-∞ω|C (ω)|2d ω=1E ∫+∞-∞ωC (ω)C 3(ω)d ω　(3表示共轭运算)52　辽宁师范大学学报(自然科学版)第32卷=1E 12π∫+∞-∞∫+∞-∞∫+∞-∞ωc 3(t )c (t ′)e j (t-t ′)ωd ωd t ′d t =1E 1j ∫+∞-∞c 3(t )d d t c (t )d t =1E ∫+∞-∞φ′(t )|c (t )|2d t(1) 从(1)可知,在整个时间范围内,对信号相位的导数关于信号时域能量密度进行积分就可以得到信号频谱的平均频率,而信号相位的导数必须是瞬时值,才可供计算平均值使用,因此将瞬时频率定义为相位的导数是很自然的.也就是说,信号频谱的平均频率等于其瞬时频率的时间平均.但对于实信号,频谱满足C (-ω)=C 3(ω),因此能量密度频谱|C (ω)2|总是关于原点对称.由于对称性,使得实信号的平均频率是0,即瞬时频率的时间平均为0,这样就无法表示出信号的物理情况.为了将瞬时频率的概念推广到实信号和多分量信号,需要解决两个问题:(1)如何将一实信号变换为一个相对应的复信号;(2)如何将一个多分量信号分解为若干个单分量信号之和.1　实信号变换为复信号的两种方法通常,我们所能获得的各种信号都是实信号,需要构造一个对应于实信号的复信号,这是由于通过复信号可以确定信号的相位,从而得到信号的瞬时频率.显然,要构造复信号可以令其实部是实信号本身,关键是如何定义其虚部,通常有两种方法:解析信号方法和正交模型方法[5].1.1　解析信号方法和正交模型方法自从Gabor 引入解析信号的概念以后,解析信号方法就成为将实信号变换为复信号的一种最常用的方法.与实信号s (t )相对应的复信号———解析信号———z a (t )定义为:z a (t )=s (t )+j H [s (t )],其中H[s (t )]表示s (t )的Hilbert 变换,即H[s (t )]=s (t )3(1/πt ).解析信号的频谱Z a (ω)在负频率部分为0,而在正频率部分是其对应的实信号的频谱在正频率部分的2倍,即:Z a (ω)=0ω<02S (ω)ω>0 由于解析信号是复信号,因此可以表示成极坐标形式,也就是用幅度和相位表示,即z a (t )=A (t )e j φ(t ).那么幅度和相位需要怎样的关系才能使得一信号是解析信号呢?幅度A (t )的频谱,S A (ω)=12π∫+∞-∞A (t )e -j ωt d t ,e j φ(t )的频谱S φ(ω)=12π∫+∞-∞e j φ(t )e -j ωt d t ,于是z a (t )=A (t )e j φ(t )的频谱Z a (ω)=S A (ω)3S φ(ω)=∫+∞-∞S A (ω-ω′)S φ(ω′)d ω′.因此,信号z a (t )的频谱Z a (ω)可看做是系数为S φ(ω′)的A (t )的已搬移频谱之和.假定S A (ω)在频率区间(-ω1,ω1)是带限的,则z a (t )是解析的,即把Z a (ω)搬移到正频率轴的充分条件是:当ω′≤ω1时,S φ(ω′)=0.因此,对于一解析信号,其低频含量在幅度上,而高频含量在e j φ(t )项上.用正交模型法将一实信号变换为复信号,首先需要将实信号写成s (t )=A (t )co s φ(t )的形式,那么其对应的复信号就是z q (t )=A (t )e j φ(t ),z q (t )称为正交模型信号.正交模型法在提出解析信号法之前就已使用.1.2　解析信号法与正交模型法的比较从表面上看,通过解析信号法得到的复信号———解析信号———似乎和正交模型法的一样都是实信号加上其正交分量,但实际上解析信号并非总是如此.这是由于在Hilbert 变换过程中首先对实信号做Fourier 变换得到双边谱,然后通过滤波得到单边谱,之后再对单边谱进行Fourier 逆变换得到复信号,因此,通过Hilbert 变换构造解析信号等价于去掉实信号频谱的负频率部分,如果实信号的正频谱泄漏到负频谱区域中,则其Hilbert 变换就不是实信号的正交分量.根据Bedrosian 乘积定理[6],如果实信号s (t )=A (t )co s φ(t )的幅值频谱Φ[A (t )]只在区间(-f 0,f 0)内有值,而频谱Φ[cos φ(t )]只在区间(-f 0,f 0)外有值,则等式A (t )cos φ(t )+j H[A (t )cos φ(t )]=A (t )e j φ(t )成立.因此,基于Hilbert 变换的解析信号发生器是一种高频选择器,信号的高频部分成为复信号的相位.它并不能准确刻画所有　第1期刘小丹等:　基于EMD的信号瞬时频率估计53实信号的物理含义,如果信号的A(t)和co sφ(t)在频域没有完全分开,则Hilbert变换就会产生一个部分重叠且相位扭曲的函数.虽然此时产生的解析信号仍是唯一的,但其结果无法预测.总之,基于Hilbert变换的解析信号法只适用于A(t)和cosφ(t)在频域完全分开的实信号,从而才能得到信号瞬时频率的较好估计.解析信号法不但有局限性,计算也比较困难,用正交模型信号来近似它可以达到相当的简化.为判断正交模型信号与解析信号在什么时候一致,需要分析这两种方法的误差,主要有两种误差度量:能量准则和逐点比较.能量准则就是通过计算实信号的解析信号z a(t)和其正交模型信号z q(t)之差的能量从整体上说明两种方法的误差,而逐点比较就是在每一时刻比较z a(t)和z q(t).无论是能量准则还是逐点比较,都可以得出结论:正交模型信号的频谱在负频率轴损失越小,则由解析信号法和正交模型法所得到的信号的一致性就越好.如果正交模型信号的频谱只位于正频率轴上,而在负频率轴上为0,则解析信号和正交模型信号是完全一致的.2　估计信号瞬时频率的正交包络法通过计算实信号的解析信号的相位导数而得到其瞬时频率,其缺点是:(1)计算解析信号比正交模型信号困难,计算量大;(2)由于解析信号通常是通过计算实信号的Fourier变换得到的,因此对于某些短促而且快速振荡的信号,会出现比较严重的G ibbs效应;(3)对于已经知道了A(t)和φ(t)的实信号,只有当A(t)和φ(t)的频谱在频域完全分隔开而没有重叠区域时,其对应的解析信号才能准确表达实信号的物理意义.对于正交模型方法,虽然其计算简单,但如何将实信号表示成s(t)=A(t)co sφ(t)的形式是一个尚未完全解决的问题[3],其实质就是怎样根据一个时间函数s(t)获得幅值A(t)和相位φ(t)这两个时间函数.在一些情况下,可以得到信号的相位或幅值,即有一个A(t)或φ(t),这样就可以使用正交模型方法来计算.但有些时候获得信号的相位或幅值非常困难.Rowe[7]提出了在使用正交模型方法表示实信号时A(t)和φ(t)必须满足两个条件:(1)s(t)= A(t)cosφ(t),A(t)≥0,该条件也可以写成s(t)=Re[A(t)e jφ(t)];(2)A(t)和φ(t)必须符合物理直觉.从条件(1)可以得出两条结论:(1)对于A(t)>0,|s(t)|=A(t)当且仅当|cosφ(t)|=1,由此可以确定A(t)和φ(t)的相切点;(2)当t从s(t)的一个零点增加到另一个零点时,φ(t)增加了π.2.1　内蕴模函数为了获得实信号s(t)的幅值A(t),我们首先对实信号s(t)本身进行限制,为此首先引入由Huang N.E.等人提出的内蕴模函数(IM F,Int rinsic Mode Function)[8]的概念.IM F是指满足以下两个条件的函数:(1)函数的过零点数目与函数的极值点数目相等或者至多相差1;(2)在任意一点,函数的上包络与其下包络的均值为0.其中,上包络是指由函数局部极大值所定义的包络,下包络是由函数局部极小值所定义的包络.Huang N.E.提出,用三次样条对IM F的局部极大值进行插值得到IM F的上包络upper;同样,用三次样条对IM F的局部极小值进行插值得到IM F的下包络lower.根据IM F的定义,显然有upper≥0且upper+lower=0.Huang N.E.等人认为,只有将信号分解成若干个IM F之和,通过分析各个IM F的瞬时频率,才能揭示原信号真正的物理意义,而且可以将这一思想应用于非平稳信号分析.2.2　经验模式分解Huang N.E.等人于1998年提出了一个自适应的、非监督的、数据驱动的多分辨分解方法:经验模式分解(EMD,Empirical Mode Decomposition)[8].该方法自适应地通过筛选过程将信号分解为局部窄带的各分量———内蕴模函数(IM F,Int rinsic Mode Function)之和,并对分解后得到的各分量IM F进行Hilbert变换,获得分量的瞬时频率和振幅,即Hilbert谱.美国NASA宇航中心将这种形式的Hil2 bert变换称为Hilbert2Huang变换,简称为H H T(Hilbert2Huang Transform)[9].EMD方法可以将非线性、非平稳过程的信号,根据信号的局部特征,自适应地分解为频率由高到54　辽宁师范大学学报(自然科学版)第32卷低的、局部窄带的各分量,即IM F.该分解算法称为筛分过程(Sifting Process ).分解模型可以表示如下:s (t )=∑pi =1c i (t )+r (t )其中s (t )为观测信号,c i (t )为第i 个IM F ,r (t )为趋势项,一般为一常值或一单调函数.该方法没有任何的能量损失,可由各分量对原信号进行重构.虽然EMD 已被广泛应用,但该方法以及使用该方法分解得到的IM F 的Hilbert 变换均存在边界效应,严重影响非平稳信号的分析.2.2.1　EMD 边界效应的抑制在EMD 方法的筛分过程中,构成上下包络的三次样条函数在数据序列的两端会出现发散现象,使边界产生较大误差,而且,这种误差随着筛分过程的不断进行而向内传播,从而“污染”整个数据序列.我们采用最大熵谱估计法,即Burg 方法来进行边界延拓.这是由于Burg 方法尤其适用于对短数据情况的预测.Burg 方法是最小化正向和反向两个预测误差的和,利用Levinson 递推关系来决定预测滤波器的参数.当用EMD 方法处理数据时,为了抑制边界效应,每次对数据进行筛分之前,我们都利用Burg 方法对数据进行延拓,具体方法是:(1)已知数据序列{x (n -p ),x (n -p +1),…,x (n -1)},用Burg 方法前向预测x (n )的值x ′(n ),产生新的数据序列{x (n -p +1),…,x (n -1),x ′(n )},同样再用Burg 方法预测x (n +1)的值x ′(n +1),依此类推,直至在前向预测值中各产生出一个新的极大值点和一个新的极小值点;(2)对同样的数据序列{x (n -p ),x (n -p +1),…,x (n -1)}用Burg 方法后向预测x (n -p -1)的值x ′(n -p -1),对新数据{x ′(n -p -1),x (n -p ),…,x (n -2)}再用Burg 方法后向预测x (n -p -2)的值x ′(n -p -2),依此类推,直至在后向预测值中各产生出一个新的极大值点和一个新的极小值点.对延拓数据得到的两个极大(小)值点和数据自身的极大(小)值点进行三次样条插值,得到上(下)包络.由此来抑制EMD 方法的边界效应,尤其是抑制EMD 方法的边界效应对低频IM F 产生的影响,这是由于EMD 方法的边界效应对高频IM F 影响较小,而对低频IM F 影响较大.2.2.2　Hilbert 的边界效应分析用EMD 方法将数据分解成IM F 和趋势项后,对各个IM F 进行Hilbert 变换,才能得到各自的瞬时频率.Hilbert 变换实际上是通过构造与原实信号具有90°相位差的共扼信号,然后构成信号的解析形式;Hilbert 变换的边界效应在求取90°共扼信号的过程中产生,共扼信号是通过“傅立叶变换—双边谱对折为单边谱—傅立叶逆变换”获得的;对于周期性信号,在非完整周期采样的情况下进行傅立叶变换,将会出现所谓的“频谱泄露”问题,将单边谱进行傅立叶逆变换的过程中,频谱泄露所造成的误差无法抵消,反映在时域波形上,将会造成所求得的共扼信号产生失真现象,这种失真主要集中在信号的两端.对于待分析的随机信号,使用Hilbert 变换法无法避免傅立叶变换过程中出现的频谱泄露问题,因此,结果存在边界效应.基于上述思想,我们提出了一个估计信号瞬时频率的新方法———正交包络法,该方法通过将IM F 的上包络来估计信号的幅值,使用正交模型法将实信号表示为复信号,未涉及傅立叶变换过程,也就没有频谱泄露的问题.2.3　正交包络法设实信号s (t )是一个IM F ,其上包络为u (t ),显然u (t )满足Rowe 提出的对信号幅值的限制,于是为了将信号表示为s (t )=A (t )cos φ(t )的形式,令A (t )=u (t ).当A (t )≠0时,令p (t )=co s φ(t )=s (t )/A (t ),对等式两边求导,有:p (t )′=-φ′(t )sin φ(t )(2) 记信号的瞬时频率为ω=φ′(t ),于是p (t )′=-ωsin φ(t ),当sin φ(t )≠0时,即|p (t )|≠1时,根据ω的非负性,有:ω=-p ′(t )sin φ(t )=|p ′(t )|1-p 2(t )(3)第1期刘小丹等:　基于EMD的信号瞬时频率估计55 当sinφ(t)=0时,即|p(t)|=1时,对式(2)两边再求导,有:p(t)″=-φ″(t)sinφ(t)-[φ′(t)]2co sφ(t)=-ω2p(t)(4)ω=p″(t)p(t)=|p″(t)|(5) 对于A(t)=0时刻的瞬时频率可以利用其他时刻的瞬时频率通过插值得到.通过以上分析,估计实信号瞬时频率的正交包络法步骤如下:1.确定信号s(t)是一个IMF;2.确定s(t)的极大值点,用三次样条对极大值点进行插值,得到s(t)的上包络u(t);3.令s(t)的幅值A(t)=u(t);4.当A(t)≠0时,令p(t)=co sφ(t)=s(t)/A(t);5.确定A(t)≠0时刻的信号s(t)的瞬时频率ω=|p′(t)|1-p2(t)当|p(t)|≠1 |p″(t)|当|p(t)|=1;6.对于A(t)=0时刻的瞬时频率可以利用其他时刻的瞬时频率通过插值得到.3　结果与讨论在实际计算中,由于信号的EMD分解存在误差,得到的IM F的上包络u(t)、下包络l(t)的均值可能不是0,为避免出现|co sφ(t)|>1的情况,修改A(t)的定义为A(t)=max(|u(t)|,|l(t)|).图1～3是Matlab的chirp函数中的3种信号,依次表示信号的瞬时频率分别是线性、二次凹函数、二次凸函数,为使实验结果比较起来更加明显,对于线性瞬时频率的chirp信号乘以一个t2.我们同时用Hilbert变换计算了信号的瞬时频率,以便两种方法的比较.需要指出的是,Matlab帮助中的chirp函数中最后一个例子有问题,按照其给出的数据无法得出其给出的结果,只有修改为如下数据才能得到相对应的结果:t=-1∶0.001∶1;f0=400,f1=100.图4中的信号是一分段信号,由两段不同频率的谐波组成,频率分别是5和10;图5中的信号是两个谐波信号与一个常数的叠加,频率分别是2和10.要估计图5中信号的瞬时频率,首先需要用EMD 算法将其分解为单分量信号,然后估计各分量的瞬时频率.　辽宁师范大学学报(自然科学版)第32卷56我们计算了使用两种不同方法得到的瞬时频率的标准差,如瞬时频率为线性函数的chirp信号,正交包络法的标准差为33.23,Hilbert变换法的标准差为36.91.从结果中可以看出,用Hilbert变换估计信号的瞬时频率存在边界效应,所得结果方差较大,而正交包络法所得到的估计边界效应很小.正交包络法的结果对信号的采样率比较依赖,在采样率较低时,无法保证采样到信号的极值点,而正交包络法正是利用信号极值点估计其幅值,于是就会产生较大的估计误差.在执行效率方面,正交包络法高于Hilbert变换法(如附表所示).在使用Hilbert变换法估计信号瞬时频率的过程中,采用解析信号法将实信号转变为复信号,这期间涉及对信号的傅立叶变换,而正交包络法直接利用在EMD分解过程中得到的上包络来估计信号的幅度,从而使用正交模型法将实信号第1期刘小丹等:　基于EMD的信号瞬时频率估计57转换成复信号,不涉及傅立叶变换,节省了运算时间,提高了执行效率.附表　使用Hilbert变换方法和正交包络法估计瞬时频率的执行时间/schirp信号的瞬时频率估计(瞬时频率为线性函数)chirp信号的瞬时频率估计(瞬时频率为二次凹函数)chirp信号的瞬时频率估计(瞬时频率为二次凸函数)分段信号的瞬时频率估计叠加信号的瞬时频率估计Hilbert变换方法0.31840.10940.23440.21881.8438正交包络法0.17190.07810.15630.09380.9531参考文献:[1]　陈平,李庆民,赵彤.瞬时频率估计算法研究进展综述[J].电测与仪表,2006,43(7):126.[2]　BOASHASH B.Estimating and interpreting t he instantaneous frequency of a signal2part1:fundamentals[J].Proc IEEE,1992,80(4):5202538.[3]　BOASHASH B.Estimating and interpreting t he instantaneous frequency of a signal2part2:algorit hms and applications[J].ProcIEEE,1992,80(4):5402568.[4]　ZHON G Y ou2ming,QIN Shu2ren.Research on t he Definition,Paradoxes and Basic Attribute of Instantaneous Frequency[C]//TAN Jiu2bin.Proceedings of t he t hird international symposium on instrumentation science and technology.Harbin:Harbin Institu2 te of Technology Press,2004:4762482.[5]　CO H EN L.Time2frequency analysis[M].Englewood Cliff s,NJ:Prentice2Hall,1995:12225.[6]　BEDROSIAN E.A product of t heorem for Hilbert transform[J].Proc IEEE,1963,51:6862689.[7]　ROWE H E.Signals and noise in communication systems[M].Princeton,NJ:Van Nostrand,1965:1082167.[8]　HUAN G N E,SH EN Z,LON G S R.The empirical mode decomposition and t he Hilbert spectrum for nonlinear and non2stationarytime series analysis[J].Proceeding of Royal Society,1998,454:9032995.[9]　ZHAO Zhi2dong,PAN Min,CH EN Yu2quan.Instantaneous Frequency Estimate for Non2stationary Signal[C]//SON G Jian.Pro2ceedings of t he5th World Congress on Intelligent Control and Automation.Hangzhou:Zhejiang University Press,2004:364123643.Estimating the instantaneous frequency of a signal based on EMDL IU Xiao2da n,　SUN Xiao2qi,　S H EN Bi n(School of Computer and Information Technology,Liaoning Normal University,Dalian116029,China) Abstract:We analyze t he concept of t he instantaneous f requency of a signal and t he two met hods of ob2 taining t he p hase of a signal,i.e.analytical met hod and ort hogonal met hod.We propo se a new met h2 od,i.e.t he ort hogo nal envelope met hod,to estimate t he instantaneous f requency based on EMD. The ort hogo nal envelope met hod is easy in calculation and it overcomes t he difficulty of ort hogonal met hod in creating two time f unctions by paring wit h Hilbert t ransform,t he ort hogonal en2 velope met hod is better in dealing wit h border effect s.The result s of experiment show t hat t he or2 t hogonal envelope met hod is effective to estimate t he instantaneous frequency.Key words:instantaneous f requency;t he ort hogonal envelope met hod;EMD;Hilbert t ransform。