基于LMD的信号瞬时频率求取方法及实验

基于LMD的信号瞬时频率求取方法及实验 DO I : 10 . 3785/ j . i ssn. 10082973 X. 2009 . 03 . 024 基于 L MD 的信号瞬时频率求取方法及实验 任达千 ,杨世锡 ,吴昭同 ,严拱标 () 浙江大学 机械与能源工程学院 ,浙江 杭州 310027 摘 要 : 研究了调频信号瞬时频率的直接求取法 ,提出了纯调频信号的瞬时频率直接求取法的适用条件 ,并用数学 方法证明了该适用条件. 针对极值点附近瞬时频率的畸变情况引进平滑处理改进了瞬时频率求取法. 应用局域均 () ( ) 值分解 L MD和经验模态分解 EMD求取仿真信号和汽轮机转子振动信号的瞬时频率. 结果 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 明 ,由 L MD 方法 求取信号瞬时频率时 ,不会出现难以解释物理意义的负频率现象. 关键词 : 时频分析 ;局域均值分解 ;瞬时频率 ;适用条件 () 文章编号 : 10082973 X 20090320523206 中图分类号 : T P206 ; T H113 . 1 ; T H165 . 3文献标识码 : A Instantaneous f requency extraction method and experiment ba sed L MD R EN Da2qia n , YA N G Shi2xi , WU Zhao2to ng , YA N Go ng2biao ( )Col l e ge o f M ec ha ni cal a n d Ene r g y En g i nee ri n g , Z he j i an g U ni ve rs i t y , H a n g z hou 310027 , Chi n a () Abstract : A new ti me2f reque ncy a nal ysi s met ho d called local mea n deco mpo sitio n L MDwa s i nt ro duce d. The re sult s of L MD a re a serial of A M2FM si gnal s , a nd t he i n st a nt a neo u s f reque ncy ca n be calculat e d by a rc co si ne f unctio n . The di rect met ho d of i n st a nt a neo u s f reque ncy e xt ractio n wa s st udied. A n app lica ble co nditio n wa s give n a nd p ro ved by mat he matical met ho d s. Co n si de ri ng t he a ber ratio n nea r t he e xt re me poi nt s , t he mo vi ng smoo t h wa s app lied to develop t he i n st a nt a neo u s f reque ncy e xt ractio n p roce ss. L MD ( ) a nd e mpi rical mo de deco mpo sitio n EMDwe re app lie d to e xt ract t he i n st a nt a neo u s f reque ncy of t he si mu2 latio n si gnal a nd vi bratio n si gnal of a t ur bi ne ro to r . The re sult s sho we d t hat t he L MD i n st a nt a neo u s f re2 que ncy wa s alwa ys po sitive Key words : time2f requency analysis ; local mean deco mpo sitio n ; instantaneo us f requency ; applicable co nditio n [ 1 ] Smit h于 2005 年提出了局域均值分解方法.信号处理 () ( local mea n deco mpo sitio n , L MD, 这种新的时频 L MD 算法和 经 验 模 态 分 解 e mpi rical mo de ) 分析方法将复杂的多分量信号分解为若干个乘积函 deco mpo sitio n , EMD方法类似 , 基于极值点来定 () 数 p ro ductio n f unctio n , P F的线性组合 , 每一 P F 义局域均值函数和局域包络函数 ,但是用滑动平均 分量由一个包络信号和一个调频信号相乘得到 ,包 代替三次样条插值. EMD 分解的最终结果是一系 络信号就是该 P F 的瞬时幅值 ,而 P F 的瞬时频率可 列本征模函数 ;而 L MD 的结果是一系列乘积函数 , 以由调频信号求出. 进一步将所有 P F 分量的瞬时 可以用反余弦函数直接求得瞬时频率. 该方法不会 频率和瞬时幅值相组合 ,即可得到原始信号的时频 出现难以解释物理意义的负频率现象. L MD 的端 [ 1 ] ( ) 点效应相比较 EMD ,在程度上轻得多 ,作用范围也分布. Smit h首先将 L MD 应用于脑电图 E E G的 收稿日期 : 2007209229 . ( ) 浙江大学学报 工学版网址 : www . jo ur nal s. zj u . edu . cn/ eng ( ) ( ) 基金项目 : 国家自然科学基金资助项目50675194;宁波市自然科学基金资助项目2007 A610014.( ) 作者简介 : 任达千19742,男 ,浙江嵊州人 ,博士生 ,从事机械故障诊断方面的研究. E2mail : boo kwo r m02 @gmail . co m 通讯联系人 :杨世锡 ,男 ,教授. E2mail : yangsx @zj u . edu . cn 比较小. 本文提出并证明了“直接法”求取信号瞬时附件 ,滑动平均的跨度取最长相邻极值点距离的 ( ) 频率的适用条件 ,与经验模态分解法进行仿真与实 1/ 3 . 设原来序列为 y i, i = 1 , 2 ,, n ,则滑动平 均的公式为验对比研究 ,表明局域均值分解有显著特点. 1 目前 ,基于 L MD 的时频分析方法虽未在机械 ( ( ) ( ) ( ) - 1+ yi= y i + N + y i + Ns 2 N + 1 故障诊断领域得到应用 ,但因一些优良的性质 :分解 ()) ) ( 3 + y i - N . 过程可避免过包络、欠包络现象 ,没有严重的端点效 式中 : 2 N + 1 即为滑动跨度 ,滑动跨度必须为奇数 ,应 , 求得的瞬时频率在任何时刻均是正的 , 基于 如果是偶数则减去 1 . 在边界附近 ,减小跨度 ,以不 L MD 的时频分析方法有望在机械故障诊断领域占 超过序列的端点为限. 一席之地. 如果在滑动平均之后 ,尚有相邻点的值相等 ,则 再次做滑动平均 , 直到任何相邻点不相等 , 如图 1 1 局域均值分解方法( ) () b和 c中的实线所示. 在计算得到局域均值函数和局域包络函数之 局域均值分解的基本思想和步骤如下 :( ) 后 ,将滑动平均后的局域均值函数 mt从原始信11 ( ) 首先确定原始信号 x t所有的局部极值点 ,分 ( ) ( ( ) ( ) ) 号 x t中分离出来 ,得到 h11 t= x t- m11 t.( ) 别由局部极值点形成局域均值函数 mt和局域11 ( ) ht除 以 滑 动 平 均 后 的 局 域 包 络 函 数再用11 ( ) ( ) 包络函数 at.原信号如图 1 a所示 ; 找出信号11 ( ) ( ) ( ) a11 t, 以 对h11 t进 行 解 调 , 得 到 s11 t= ,任意 2 个相邻极值点之平均值即定的极值点之后 ( ) ( ) ( ) 如果 s11 t相对应的局域包络函数ht/ at.11 11 ( ) 义为局域均值 ,如图 1 b虚线所示.第 i 段的局域( ) ( ) a12 t= 1 , 则已经是一个调频信号. 如果 a12 t? 均值函数由下式给出 :( ) 1 ,则将 st作为原始数据重复以上迭代过程 ,直11 ( ) ()1 m = ni + ni +1 / 2 .i( ) ( ) 到 s1 n t?1 ,它的1 n t 为一个调频信号 ,即 - 1 ?s 式中 : ni 和 n i +1 为相邻极值点.( ) 局域包络函数 at = 1 . 最后把迭代过程中产( ) 1 n+1 同样 ,第 i 段局域包络函数为相邻两极值点差 生的所有局域包络函数相乘便可以得到 P F 的包络 n () 值的 1/ 2 ,如图 1 c虚线所示. ( ) 信号 at( ) ( ) ( ) ( ) ,该1 atatatat= 11 12 1 n = 1 q 7 第 i 段局域包络函数由下式给出 :q = 1 ()a= n- n/ 2 . 2 ( ) i i i +1 包络信号即 P F 分量的瞬时幅值. 将包络信号 a1 t ( ) 和调频信号 st相乘便可得到原始信号的第一个然后在此基础上 ,用同样的算法对局域均值函1 n ( ) ( ) ( ) P F 分量 Y t= atst,它包含了原始信号 PF1 1 1 n 数和局域包络函数分别作滑动平均.据文献 [ 1 ] 的中最高的频率成分 , 是一个单分量的调频2调幅信 ( ) 号 ,其瞬时幅值就是包络信号 at, 其瞬时频率 1 ( ) ( ( ) ) f t则可由调 频信 号 st求出. 将 x t减 去 1 1 n ( ) ( ) ( ) Yt即得到一个新的信号 ut,将 ut作为原 PF1 1 1 ( ) 始数据重复以上步骤 ,循环 k 次 ,直到 ut为一个单 1 ( ) 调函数为止 ,这样便可以将 x t分解为 k 个 PF 分量 和一个单调函数 u之和 ,即k k ( ) ( ) ( ) x t= Y t+ ut.PF n k ?n = 1 2 “直接法”求取调频的瞬时频率 L MD 分解得到的 P F 是一单分量的调频调幅 ( ) 信号 ,其瞬时幅值就是包络信号 a t,其瞬时频率 φ( ) ( ) ( ) f t则可由调频信号 s t= co s t求出 ,即 ()φ( ) ( ( ) ) 4 t= a rcco s s t. 图 1 局域均值函数和局域包络函数( ) ( ) φ 将t展开然后求导 ,即可求得 s t的瞬时频Fig. 1 Local mean f unctio n a nd local envelop f unctio n 率 ,也就是 P F 的瞬时频率 ,即 ()() () ω( ) φ( ) 5 t= dt/ d t. 将式 4代入式 7,可得 ( ) - 1 特别需要指出的是 ,必须保证 - 1 ?s t?1 , ()( ) ( ) 8 f t= s′t.2 ( ) π ( ) 21 - sts t的极点应等于 ?1 ,如果不等于而接近于 ?1 , ( ) ( ) 当 s t= ?1 时 ,式 8的分母等于 0 . 因此这则可设置成 ?1 . ( ) ( ) 是一个不定式 , s′t的误差将被放大. 在 s t= ?上述求取瞬时频率的方法显得简单和直接 ,称 ( ) 1 的一个邻域内 , 式 8并不能得到可靠的计算结之为“直接法”. “直接法”的一个优点是计算得到的 () 果 ,这就是图 3 c所示的瞬时频率计算结果出现明 瞬时频率都是正的 ,不会出现难以解释物理意义的 显毛刺的原因. 因此 ,必须作一些处理 ,将这一些点 负频率. 去除. 以邻近点的平均值代替或作平滑处理 ,得到下面结合实例来说明“直接法”求取瞬时频率的 的瞬时频率如图 4 所示.计算过程. 例 1 调频信号 ,波形如图 2 所示 ,幅值无物理 单位 ,表达式为 2 ( ) () () s t= co s 50 t + 40 t; 0 < t ?1 .6 π() 因反余弦函数的值域为 [ 0 ,] ,将式 6代入式 () ( ) 4计算得到的相位如图 3 a 所示 ; 相位展开成可 图 4 平滑处理后瞬时频率( ) 导的单调递增函数 ,如图 3 b所示. Fig. 4 Smoo t hed instanta neo us f requency 相位展开后再对求其导数 , 可得到瞬时频率. ω( ) ( π) 为方便后续处理 , t乘以 1/ 2,得 3 “直接法”适用条件φ( ) 1 dt 1 ( ) ω( ) ()t= t=f ×.7 π π 22d t () , 原始信号经 L MD 分解“直接法”求瞬时频率 图 3 c即为计算所得的瞬时频率. 可以看到图中 ( ) ( ) 有一些明显的毛刺 ,出现的位置是在 s t= ?1 附近.得到的调 频信 号 s t应满 足一 定的 适 用 条 件. [ 1 ] Smit h提出信号必须在[ - 1 ,1 ] 的范围内 ,这是因 为 a rcco s 函数的定义域即为[ - 1 ,1 ] . 从瞬时频率 ( ) 的概念出发 , 为得到正的瞬时频率 , 调频信号 s t 必须满足更苛刻的条件. 应用“直接法”求取信号的 瞬时频率 ,信号必须满足如下条件 : 信号的极大值等于 1 ,信号的极小值等于 - 1 . ( ) φ( ) ω证明 : 为保证瞬时频率 t= dt/ d t 有定 φ( ) ω( ) φ( )义 ,t必须是可导的 ,为了使t> 0 ,相位t 必须是连续单调递增函数. [ 3 ] φ( ) 根据 dt/ d t 的符号和三角函数的性质: φ( ) (φ( ) π) co s t= co s t+ 2 m, ()9 φ( ) ( φ( ) π) co s t= co s - t+ 2 n. 式中 : m 、n 为任意整数. ( ) φ 以t的极值点为界将原相位函数分段 ,重新 定义一相位函数 ,得 φ( ) πφd/ d t > 0 ; t+ 2 n, Φ( ) φ( ) πφ() d/ d t 10 t= - t+ 2 m,< 0 ; 其他.未定义 , ( ) Φ 如果 t可以构成一连续单调递增函数 , 即 φ( ) ( ) ( ) Φφt的极值点是 t的可去间断点 ,那么t即 ( ) 是可以展开的 ,如例 1 的图 3 b. Φ( ) t的连续性 , 即由间断点的特性决定. 在 φ( ) φ( ) t的极大值点处 , 左导数 ′t> 0 , 右导数- φ( ) ( ) Φ′t< 0 ,那么 t的左右极限为+ - ( ) Φ( ) φ( ) π6 b所示. t= t+ 2 n,()11 + Φ( ) φ( ) π若以此仿真信号作为原始信号 , 经 L MD 分解 t= - t+ 2 m. ( ) ( ) () 结果如图 7 a 、b和 c所示 ,图中信号的幅值无 Φ( ) φ( ) t在t的极大值点处连续的充要条件为- + Φ( ) Φ( ) t= t,因此有 物理单位. 再进一步用“直接法”可求得瞬时频率 , φ( ) t= 0 , ( ) 如图 7 d所示. ()12 n = m . 或 φ( ) π t= ,()13 n = m - 1 . 为此信号的极大值点可用下式表示 : ()( ) (φ( ) ) 14 s t= co s t= ?1 . ( ) φ 当t为极小值点时 ,可以得到同样的结论 , 即信号的极值点必须是 ?1 . “直接法”的适用条件 ,有如下意义 : ) 1给出了信号应用“直接法”求取瞬时频率必 图 6 “直接法”计算所得相位须满足的条件 ,为进一步研究调频信号的判据提供 Fig. 6 Pha se generated by direct met ho d 了理论基础. ) 2说明瞬时频率不仅仅是局部的概念 ,为了计 5 与经验模态分解的实验对比算某一时刻的瞬时频率 ,至少需要知道信号的全部 5 . 1 经验模态分解简介 极值点. [ 2 ] H ua ng 等 人提 出 了 一 种 被 称 为 Hil ber t2 ( ) H ua ng 变换 简称 H H T的信号分析方法 , 其核心 4 “直接法”的仿真研究 ( 内容是经验模态分解 e mpi rical mo de deco mpo si2 () tio n , EMD,将信号分解成若干个本征模函数 i n2 以下例子可以说明”直接法”的适用条件. ) t ri n sic mo de f unctio n , IM F. 对各 IM F 做 Hil be r t 例 2 仿真信号波形如图 5 所示 , 幅值无物理 变换 ,可以得到信号的瞬时频率 ,并进一步计算信号 单位 ,信号表达式为 2的时频分布和边缘谱. ( ) ( ) ()x t= co s at + bt + c. 15 为了得到有物理意义的瞬时频率 , H ua ng 等 式中 : a = 1 , b = - 10 , c = 10 .[ 6 ] 人提出了归一化 Hil ber t 变换. 从 EMD 分解得 在 t = 5 处 ,信号有一极小点 ,其值为 - 01759 7. 到的 IM F 出发 ,首先找出 IMF 的所有极大值点 ,对 信号的相位表达式为( ) 这些极大值点样条插值得到包络线 Et;再将 IMF 2 φ( ) ()t= at + bt + c. 16 ( ) 除以 Et即得到归一化的 IM F ; 然后对归一化的 信号的瞬时频率表达式为IMF 进行 Hilbert 变换求瞬时频率. Huang 认为归一 ω()( ) φ( ) t= dt/ d t = 2 at + b. 17 化 Hilbert 变换的结果是令人满意的 , 但并非所有 但是“直接法”不能计算此信号的瞬时频率 ,按 IMF 都可以求得有物理意义的瞬时频率. () ( ) 式 4计算得到的相位如图 6 a所示. 因为信号在 5 . 2 仿真信号对比t = 5 处有一极小值点不等于 - 1 ,不满足“直接法” L MD 方法的分解过程可以得到一系列的 P F( ) Φ的适用条件 ,相位连续展开后 t在 t = 5 处有一 分量 , P F 分量由包络信号和调频信号相乘得到 ,从 间断点 ,所以不能得到有物理意义的瞬时频率 ,如图调频信号计算得到的瞬时频率是正的、连续的、具有 物理意义的. 而 EMD 方法先得到 IM F 分量 ,再对 IM F 分量进行 Hil be r t 变换求得瞬时频率和瞬时幅 值 ,则可能产生无法解释的负频率. 以图 8 所示信号为例 ,信号的幅值无物理单位 , 图 5 仿真信号表达式为 Fig. 5 Simulated signal 图 8仿真信号 Fig. 8 Simulated signal 的过渡处出现波动 ,还出现负的瞬时频率值 ,并无法 解释其物理意义 ,如图 9 所示. 图 10 所示为应用“直接法”求取的瞬时频率 ,与 Hil ber t 变换最大的不同是瞬时频率总是正的 ,没有 出现无法解释的负频率. [ 425 ] 负频率现象不仅是一个理论问题,而且在分 析实际信号时也是经常出现的 ,成为瞬时频率估计 的一个重要的误差来源. 5 . 3 振动信号对比 如图 11 所示 ,这是一段从汽轮机上采集到的转 子振动信号 , EMD 分解的结果见图 12 ,信号幅值的 物理单位是 mV . 可以看到总共有 9 个 IM F 和一个 趋势项 ,因端点效应的影响 IM F9 和趋势项的幅值 远超出了原信号 ,部分信号无法显示. 图 13 显示了 IM F1 的瞬时频率 ,可以清楚地看到负频率现象. 图 14 为振动信号的 L MD 分解结果 ,信号幅值 图 7 仿真信号 LMD 分解结果及瞬时频率 Fig. 7 L MD result f ro m simulated signal a nd in2 sta nta neo us f requency ( π ) co s 2×5 t,0 ?t < 0 . 63 ; ( ) x t= - 1 . 02 x + 1 . 235 4 , 0 . 63 ?t < 0 . 635 7 ; ( π ) co s 2×5 t - 0 . 2,0 . 635 7 ?t ?1 . ()18 信号前后两段均为余弦信号 ,在 01 63 s 附近信号 的相位有一个突变 ,前后信号相位差 01 2 rad ,在连接 处 ,用一段直线过渡. 这个函数既符合 IMF 的定义 , 可以用 Hilbert 变换求取瞬时频率 ,也符合“直接法” 的适用条件 ,可以用“直接法”求取瞬时频率. 该信号经 Hil ber t 变换得到的瞬时频率在函数 了理论基础. () 2针对极值点附近畸变情况引进平滑处理改 进了该算法. () 3与经验模态分解法进行了实验对比 ,表明局域 均值分解法具有避免出现负频率现象的显著特点. ( ) 参考文献 Ref erences: [ 1 ] SM I T H J S. The local mea n deco mpo sitio n a nd it s ap2 plicatio n to E E G percep tio n data [ J ] . Journal of The () Royal Society Interface , 2005 , 2 5: 4432454 . [ 2 ] H U A N G N E , S H EN Z , L ON G S R , et al . The empir2 ical mo de deco mpo sitio n a nd t he Hil bert spect r um fo r 的物理单位是 mV . 图 15 是用“直接法”求取 P F1 no nlinear a nd no n2statio na r y time series analysi s [ J ] . 的瞬时频率 ,可以看到用“直接法”求得瞬时频率总 ( ) Proceedings of the Royal Society A , 1998 , 454 1971 : 是正的 ,不会出现难以解释物理意义的负频率现象. 9032995 . [ 3 ] 中国矿业大学数学教研室. 数学手册[ M ] . 北京 : 科学 出版社 , 1990 : 73295 . 6 结论[ 4 ] 科恩 L . 时2频分析 : 理论与应用 [ M ] . 白居宪 , 译. 西 安 : 西安交通大学出版社 , 1998 : 22233 . 局域均值分解是一种新出现的时频分析方法 ,[ 5 ] BOA S HA S H B . Estimating a nd i nterp reting t he insta n2 所求得的瞬时频率不会出现难以解释物理意义的负 ta neo us f requency of a signal2pa rt 1 : f unda mental s [ J ] . () 频率. 本文的分析研究有以下几方面的意义 : Proceedings of the IEEE , 1992 , 80 4: 5202538 . H U A N G N E , A T TO H2O KIN E N O . The Hilbert2 [ 6 ] () 1给出了信号应用“直接法”求取瞬时频率必 Huang transf orm in engineering [ M ] . Boca Rato n : 须满足的条件 ,为进一步研究调频信号的判据提供 Taylo r &Fra nci s Gro up , 2005 .