用希尔伯特黄变换(HHT)求时频谱和边际谱

Python希尔伯特黄变换（Python Hilbert-Huang Transform，简称HHT）是一种复杂非线性信号分析方法，结合了希尔伯特变换和黄变换的优势，能够有效地对非线性和非平稳信号进行时频谱分析。

本文将从HHT的原理、基本步骤和Python实现方法三个方面进行介绍。

一、HHT的原理1.希尔伯特变换希尔伯特变换是一种将实数信号转换为解析信号的数学方法，通过对原信号进行傅立叶变换得到频谱信息，再对频谱信息进行一定的处理得到解析频谱，从而实现信号的解析表示。

希尔伯特变换的核心是求出原信号的解析函数，即原信号的复数形式，其中实部是原信号本身，虚部是原信号的希尔伯特变换。

希尔伯特变换在信号处理领域有着广泛的应用，能够提取信号的瞬时特征，对非平稳信号进行时频分析具有很高的效果。

2.黄变换黄变换是一种局部线性和非线性信号分解方法，可以将非线性和非平稳信号分解成若干个固有模态函数（Intrinsic Mode Function，简称IMF）的线性组合。

黄变换首先对原信号进行极值点的提取，然后通过极值点之间的插值得到包络线，再将原信号减去包络线得到一维信号，并对得到的一维信号进行数据挑选和插值，最终得到IMF。

多次重复以上步骤，直到原信号能够被分解为若干个IMF，再通过IMF的线性组合得到原信号的近似表示。

3.HHT的结合HHT将希尔伯特变换和黄变换结合在一起，利用希尔伯特变换提取信号的瞬时特征，再通过黄变换将信号分解成若干个IMF，从而能够更准确地描述信号的时频特性。

HHT的优势在于能够适用于非线性和非平稳信号，对信号的局部特征具有很好的描述能力，因此在振动信号分析、生物医学信号处理等领域有着广泛的应用。

二、HHT的基本步骤1.信号分解HHT首先对原信号进行希尔伯特变换，得到信号的瞬时频率特征，然后通过黄变换将信号分解成若干个IMF。

2.IMF的提取针对得到的IMF，需要对每个IMF进行较为严格的判别，确定其是否符合IMF的特征：极值点交替出现、包络线对称、局部频率单调。

希尔伯特黄变换及其应用希尔伯特黄变换及其应用希尔伯特黄变换（Hilbert-Huang Transform，HHT）是一种用于分析非线性和非平稳信号的方法，它由黄其森（Norden E. Huang）和希尔伯特（Hilbert）共同提出。

该方法通过将信号分解为一组固有模态函数（Intrinsic Mode Functions，IMF）来提取信号中的模式和趋势。

本文将介绍希尔伯特黄变换的应用，并详细讲解其中的几个应用领域。

应用一：信号处理•希尔伯特黄变换可以用于音频信号处理，通过提取信号的固有模态函数，可以分离出音频信号中的主要频率成分，从而实现去噪、降噪等处理。

•在图像处理中，希尔伯特黄变换可以用于边缘检测和纹理分析。

通过提取图像的固有模态函数，可以分离出图像中的纹理信息和边缘信息，从而实现图像增强和分割等操作。

应用二：地震学•地震学中的信号分析是一项重要的任务，希尔伯特黄变换可以用于地震信号的分析和处理。

通过将地震信号分解为固有模态函数，可以提取出地震信号中的地震波的时频特征，从而实现地震信号的分类和识别。

•希尔伯特黄变换还可以用于地震信号的时频谱分析，通过将地震信号分解为固有模态函数，并对每个分量进行傅里叶变换，可以得到地震信号的时频谱图，从而更好地理解地震信号的时频特性。

应用三：医学工程•在医学工程中，希尔伯特黄变换可以用于生物信号的分析和处理，如心电图（ECG）和脑电图（EEG）等。

通过将生物信号分解为固有模态函数，可以提取出信号中的重要特征，如心跳频率、脑电波的频率等，从而实现疾病的诊断和监测。

•希尔伯特黄变换还可以用于生物信号的时频谱分析，通过将生物信号分解为固有模态函数，并对每个分量进行傅里叶变换，可以得到信号的时频谱图，从而更好地分析信号的时频特性。

应用四：金融市场•在金融市场中，希尔伯特黄变换可以用于股票价格的分析和预测。

通过将股票价格分解为固有模态函数，可以提取出股票价格的趋势和周期成分，从而更好地预测股票价格的走势。

希尔伯特-黄变换（Hilbert-Huang Transform，HHT）0 前言传统的数据分析方法都是基于线性和平稳信号的假设，然而对实际系统，无论是自然的还是人为建立的，数据最有可能是非线性、非平稳的。

希尔伯特-黄变换（Hilbert-Huang Transform，HHT）是一种经验数据分析方法，其扩展是自适应性的，所以它可以描述非线性、非平稳过程数据的物理意义。

1 HHT简介[贺礼平.希尔伯特-黄变换在电力谐波分析中的应用研究[D].湖南:中南大学,2009]HHT的发展。

1995年，Norden E.Huang为研究水表面波构思出一种所谓“EMD--HSA”的时间序列分析法，通过这种方法他发现水波的演化不是连续的，而是突变、离散、局部的。

1998年，Norden E.Huang等人提出了经验模态分解方法，并引入了Hilbert谱的概念和Hilbert谱分析的方法，美国国家航空和宇航局（NASA）将这一方法命名为Hilbert-Huang Transform，简称HHT，即希尔伯特-黄变换。

HHT是一种新的分析非线性非平稳信号的时频分析方法，由两部分组成：第一部分为经验模态分解（Empirical Mode Decomposition，EMD）（the sifting process，筛选过程），它是由Huang提出的,基于一个假设：任何复杂信号都可以分解为有限数目且具有一定物理定义的固有模态函数（Intrinsic Mode Function，IMF；也称作本征模态函数）；EMD方法能根据信号的特点，自适应地将信号分解成从高到低不同频率的一系列IMF；该方法直接从信号本身获取基函数，因此具有自适应性，同时也存在计算量大和模态混叠的缺点。

第二部分为Hilbert谱分析（Hilbert Spectrum Analysis，HSA），利用Hilbert变换求解每一阶IMF 的瞬时频率，从而得到信号的时频表示，即Hilbert谱。

收稿日期:2002-05-30基金项目:国家自然科学基金资助项目(50178055)作者简介:罗奇峰(1948-),男,江西南昌人,研究员,工学博士,博士生导师.E -mai l:luo@mail.tongji Hilbert -Huang 变换理论及其计算中的问题罗奇峰,石春香(同济大学结构工程与防灾研究所,上海 200092)摘要:以地震波的谱分析为例,对比分析了Hilbert -Huang 变换(HHT )与傅立叶变换的差别,结果表明:HHT 是对非平稳时程进行数据分析的有效工具,能有效地将各种频率成份以固有模态函数的形式从时程曲线中分离出来,但这种固有模态函数与傅立叶变换结果不同,而且Hilber t 谱是包含时间-频率-振幅的三维谱;可采用在端点改造一个小波串的方法解决HHT 存在的端点飞翼现象.文中还提出确定是否终止固有模态函数分离过程的一种方法.关键词:Hilbert -Huang 变换;固有模态函数;傅立叶变换;端点飞翼中图分类号:T U 352 文献标识码:A 文章编号:0253-374X(2003)06-0637-04Hilber-t Huang Transform and Several Problem s inIts Calculation MethodL UO Qi -f eng,SHI Chun -x iang(Research Institute of Stru ctural Engineering and Di saster Reducti on,Tongji University,Shanghai 200092,China)Abstract :H ilbert -Huang Transform (HH T)method,w hich includes H uang transform and H ilbert transform,is applied to make spectra analysis of seism ic waves.The comparison of HHT and Fourier Transform shows that HHT is an efficient method for analysis of non -stationary data.In HHT different intrinsic mode functions (IMF)can be separated by using empirical mode decomposition (EM D),i.e.H uang Transform,from time history.IM F is different from Fourier spectrum,and H ilbert spectrum is a tempora-l frequency -amplitude spec -trum which includes three dimensions.The HHT process show s that the suggested method,by improving one short wave series at the end of the time history ,can solve the end sw ing problem.A method to end EM D pro -cess is also suggested in this paper.Key words :Hilbert -Huang transform;intrinsic mode function;Fourier transform;end sw ing近年来发展起来的H ilbert -H uang 变换(HHT )理论,是一种适合分析非平稳时程的谱分析方法[1].本文首先介绍HHT,然后以地震波的谱分析为例,对比分析HH T 和傅立叶谱分析方法,并对在实现HHT 中出现的问题进行分析与讨论.1 Hilbert -Huang 变换H ilbert -Huang 变换由Huang 变换和H ilbert 谱分析两部分组成.1.1 H uang 变换Huang 变换的关键是经验模态分离法.该方法认为任何复杂的时间序列都是由一些相互不同的、简单第31卷第6期2003年6月同 济 大 学 学 报JOURNAL OF T ONGJI UN IVERSIT Y Vol.31No.6 Jun.2003的、并非正弦函数的固有模态函数组成.基于此可从复杂的时间序列直接分离成从高频到低频的若干阶固有模态函数,即基本时间序列.固有模态函数需满足以下两个条件[1]:¹对整个时间序列来说,极值的个数与穿过零点的个数相同或其差值为1;º在任何一点最大值包线与最小值包线的均值为零.图1为一条原始的地震记录.图2为经Huang变换后得到的一个典型的固有模态函数,从图中曲线可看出其极值的个数与曲线穿过零点的个数相同,且该时程的上下包线相对零线对称.如图2所示,对任意一条时程曲线都可按下述方法进行H uang变换[1,2].图1一条原始的地震记录Fig.1A seismicrecord 图2经Huang变换后得到的一个典型的固有模态函数Fig.2A typical intrinsic mode function首先,确定时程曲线X(t)的所有峰值点,然后用三次样条函数曲线循序连接所有的最大值,得到时程曲线X(t)的上包络线.采用同样的方法循序连接所有的最小值,得到时程曲线X(t)的下包络线.循序连接上、下包络线的均值可得一条均值线m1(t),于是可得h1(t)=X(t)-m1(t).如果h1(t)满足固有模态函数所需的条件,则h1(t)即为第一阶固有模态函数.一般地说,它并不满足条件,此时,将h1(t)看成新的时间曲线,重复上述方法,可得h11(t)=h1(t)-m11(t).这里,m11(t)是h1(t)的上、下包络线.按上述方法重复k次,直到h1k(t)满足固有模态函数的条件为止.h1k(t)由下式计算:h1k(t)= h1(k-1)(t)-m1(k-1)(t).式中的h1k(t)即为第一阶固有模态函数,记作c1(t),c1(t)=h1k(t).第一阶固有模态函数c1(t)包含着时程X(t)的频率最高的成份.从X(t)中减去高频成分c1(t),得到频率较低的残差为r1(t)=X(t)-c1(t).将r1(t)看成一组新的数据,重复上述经验模态分离过程.经过多次运算可以按下式得到所有的r j(t),r j(t)=r j-1(t)-c j(t),j=2,3,,,n.当满足以下两个条件之一时,整个振型分离过程可以终止.这两个条件是:¹c n(t)或r n(t)小于预定的误差;º残差r n(t)成为一单调函数,此时不可能再从中提取固有模态函数.最后时程曲线X(t)可以按式(1)表示成n阶固有模态函数和第n 阶残差r n(t)之和.X(t)=E n j=1c j(t)+r n(t)(1)图3为对图1中的地震波作Huang变换得到的不同阶的固有模态函数.对比傅立叶变换不难看出:傅立叶变换得到的是一系列具有相同振幅的、单一频率的时程曲线;而通过Huang变换得到的是一系列固有模态函数,每一固有模态函数含有不同频率成分,且每一时刻的振幅也不尽相同,固有模态函数的阶数愈低、其所含高频成份愈多.因此可以说H uang变换比傅立叶变换更具一般意义.1.2H ilbert谱分析对给定的连续时程曲线X(t),其Hilbert变换定义为[3]Y(t)=1P Q X(S)t-S d S(2)则有Z(t)=X(t)+i Y(t).Z(t)称为X(t)的解析信号,可写为Z(t)=a(t)ex p[t H(t)](3)其中:a(t)=X2(t)+Y2(t),H(t)=arctan[Y(t)/X(t)].按式(3)的极坐标表达式,瞬时频率可定义为X(t)=d H(t)/d t[1].可见,瞬时频率X(t)也是时间的函数,对X(t)的n阶固有模态函数c(t)进行H ilbert变换,则H(t)为638同济大学学报第31卷图3 Huang 变换得到的固有模态函数Fig.3 Intrinsic mode f unctions from Huang transf ormH (t)=E nj=1a j (t)exp i Q X j (t )d t (4)其中a j (t)是第j 阶固有模态函数c j (t )的解析信号的幅值.对照式(1),这里省略了第n 阶残差r n (t),这是因为r n (t )是单调函数或常数的缘故.H (t )也可以表示为以下傅里叶变换形式[1]:H (t)=E n j=1a j exp (i X j t)(5)式中的a j 和X j 都是常数.对比式(4)和式(5),两者有相似的形式,可以说式(4)是一种广义的傅里叶变换.实际上,式(4)的H (t)既是时间的函数,又是瞬时频率X 的函数,而瞬时频率X 也是时间的函数,因此,取式(4)的实部,定义它为Hilbert 谱,记作H (X ,t).图4a 直接将H ilbert 谱绘在时间)频率)振幅的三维坐标系中,图4b 是用等值线将Hilbert 谱表示在时间)频率的两维坐标系中.从地震波的H ilbert 谱中可以看出振幅即能量随时间和瞬时频率分布的特征,与时间)频率反应谱有相同作用[4],这对研究地震波的工程特征具有重要意义.图4 Hilbert 谱的时间)频率)振幅三维图Fig.4 Hilbert spectra in tempora-l frequency -amplitude three dimensional coordinates639 第6期罗奇峰,等:Hilbert -Huang 变换理论及其计算中的问题2 HHT 变换中出现的问题及解决方法2.1 固有模态函数分离标准的选择如前所述,Huang 变换中的经验模态函数分离法的本质是筛选,即从数据中分离出满足条件的固有模态函数.满足固有模态函数的第一个条件,可以消除附加波的影响;而满足第二个条件常常是难以做到的,需要确定一个标准使得这一分离过程能够停下来.H uang 等提出通过限制标准差S [1]的大小来确定,即:S =E n k =1h 1(k-1)(t)-h 1k (t)2h 21(k-1)(t).其中:S 值定在0.2与0.3之间.本文是用分离结果的最大值的包线与最小值的包线的均值是否小于给定的小数值,来确定是否终止固有模态函数分离过程.固有模态函数分离终止标准取得不同,分离出的固有模态函数的个数和振幅也各异.2.2 H HT 中的端点飞翼问题在HHT 中会出现两种端点问题,一种是出现在三次样条拟合中,另一种是出现在H ilbert 变换中.在样条拟合中最严重的问题就发生在端点上,如果在端点处不对三次样条进行处理,就会在固有模态函数分离过程中出现很大或很小幅值(见图5)的现象,称之为端点飞翼.对此问题,可以通过使用一些改进的样图5 三次样条端点飞翼问题Fig.5 End swing phenomena at the ends of cubic splines 条函数方法,如收紧样条函数[5,6]加以解决.本文采用的方法在两个端点各改造一个波串,首先寻找数据端部的两个极值点,确定极值间的点数;用最后两个连续极值间(极大值和极小值之间或极小值和极大值之间)的波的1.5个周期成份组成一个波串,用此波串替代原数据端部相同数据长度的波.用这种方法可以较好地消除在端点出现的幅值/飞翼0现象.3 结论本文阐述了HHT 理论和它与傅立叶谱的区别,对实现HHT 的过程中存在的问题进行了讨论与分析,提出了解决办法,得出以下结论:(1)本文的工作也说明HHT 能将各种频率成份以固有模态函数的形式,从地震波的非平稳时程信号中分离出来,但这种固有模态函数与傅立叶谱不同,每一阶固有模态函数都不是单色的;Hilbert 谱是包含时间)频率)振幅的三维谱.(2)固有模态函数分离终止标准不同,分离出的固有模态函数也有差异.本文是用分离结果的最大值的包线与最小值的包线的均值是否小于一给定的小数值,来确定是否终止固有模态函数分离过程.(3)HHT 存在的端点飞翼现象,可以采用在端点改造一个小波串的方法解决.本文对HHT 的研究是初步的,它与小波理论的对比分析等问题将在另文讨论.致谢:美国科罗拉多矿业大学的张瑞冲博士、马硕硕士提供了宝贵的资料与帮助,在此表示衷心的感谢.参考文献:[1] Huang N E,Shen Z,Long S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non -station ti m e seriesanalysis[J].Proc R Soc,1998,A454:903-995.[2] S huo M.HHT analysi s of near -field seismic ground motions[D].Denver:Colorado S chool of M ines,2001.[3] 胡广书.数字信号处理)))理论、算法与实现[M ].北京:清华大学出版社,1997.[4] 张晓哲,罗奇峰.地震加速度波的时-频反应谱分析[J ].结构工程师,1999,(增刊):37-42.[5] de Boor C.A practical guide to splines[M ].New York:Springer Verlag,1978.[6] S chumaker L L.S pline functions:Basic theory[M ].New York:W i ley,1981.640 同 济 大 学 学 报第31卷。

用希尔伯特黄变换（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('瞬时频率')小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。

希尔伯特黄变换算例2电⼒⼯程信号处理应⽤希尔伯特黄变换【⽬的】1．了解希尔伯特黄变换的理论知识及应⽤领域2．⽤Matlab软件仿真，验证希尔伯特黄变换的优点【希尔伯特黄变换】希尔伯特黄变换（Hilbert-Huang transform, HHT）⾸先采⽤EMD⽅法将信号分解为若⼲个IMF分量之和，然后对每个IMF分量进⾏Hilbert变换得到的瞬时频率和瞬时幅值，从⽽得到信号的Hilbert谱，Hilbert谱表⽰了信号完整的时间-频率分布，是具有⼀定的⾃适应的时频分析⽅法。

与前⾯的⼩波分析⽅法相⽐，避免了⼩波分析基选取的困难。

分析⾮线性、⾮平稳信号采⽤基于经验模态分解的HHT⽅法可以较好地分析信号的局域动态⾏为和特征。

由于HHT⽅法的种种特点，其在机械振动、⽣物医学、故障诊断、海洋学科、地震⼯程学以及经济学各学科中得到了⼴泛应⽤。

在电⼒系统领域中，HHT⽅法可⽤于谐波分析、同步电机参数辨识、低频震荡分析、电能质量检测、磁铁谐振过电压辨识等⽅⾯和超⾼速⽅向保护等⽅⾯。

HHT⽅法在电⼒系统中的应⽤还在进⼀步的研究和探索中。

【EMD 分解】对于⼀个时间序列()x t ,其经验模态分解过程如下：（1）确定原始信号()x t 的所有极⼤值点和极⼩值点；（2）采⽤样条函数求出()x t 的上、下包络线，并计算均值()m t ；（3）做差()()()h t x t m t =-；（4） ()h t 是否满⾜终⽌条件，若不满⾜将()h t 作为新的输⼊信号转⾄第（1）步，否则转为第（5）步；（5）令()c h t =，c 即为⼀个IMF 分量，做差()r x t c =-；（6） r 是否满⾜终⽌条件，若不满⾜则将r 作为新的输⼊信号转⾄第（1）步，若满⾜则EMD 分解过程结束，不能提取的为残余量。

具体流程如图1所⽰。

EMD 分解过程图1 EMD 分解流程图对于分解总阶数为n 的时间序列，最后可以表⽰成1()()()ni i x t c t r t ==+∑式中，()r t 为残余函数，它是以单调函数。

黄锷院士在《On Holo-Hilbert spectral analysis: a full informational spectral representation for nonlinear and non-stationary data》中提出一种高维全息谱分析理论HHSA(Holo-Hilbert spectral analysis)要理解HHSA方法，首先要了解希尔伯特变换、经验模态分解(EMD)、与希尔伯特-黄变换(HHT)。

学术背景：在信号处理与频谱分析的目的是要描述信号的频谱含量在时间上变化，以便能在时间和频谱上同时表示信号的能量或者强度。

傅里叶频谱并没有告诉我们哪些频率在什么时候出现。

因此傅里叶变换无法表现信号频率成分的时变性，因此学术界先后发展出了短时傅里叶变换、窗口傅里叶变换、小波等手段，近似的求信号某一时刻的瞬时频率。

希尔伯特变换：希尔伯特变换是以著名数学家大卫·希尔伯特(David Hilbert)来命名。

通过希尔伯特变换，使得我们对短信号和复杂信号的瞬时参数的定义及计算成为可能，能够实现真正意义上的瞬时频率的提取，因而希尔伯特变换在信号处理上具有十分重要的地位，使得希尔伯特变换具有广泛的工程应用。

但在进一步的工程应用中，希尔伯特变换具有以下缺陷：(1)希尔伯特变换只能近似应用于窄带信号。

但实际应用中，存在许多非窄带信号，希尔伯特变换对这些信号无能为力。

即便是窄带信号，如果不能完全满足希尔伯特变换条件，也会使结果发生错误。

而实际信号中由于噪声的存在，会使很多原来满足希尔伯特变换条件的信号无法完全满足；(2)对于任意给定时刻，通过希尔伯特变换运算后的结果只能在一个频率值，即只能处理任何时刻为单一频率的信号；(3)对于一个非平稳的数据序列，希尔伯特变换得到的结果很大程度上失去了原有的物理意义。

图1 傅立叶、小波与希尔伯特-黄变换对瞬时频率的分辨率希尔伯特-黄变换：针对上述的三个问题，黄锷院士在1998年提出希尔伯特-黄变换(HHT)。