基于R实现Hilbert—Huang变换的网络技术分析

　1 简介
　　数字信号处理技术应用广泛，包括金融、通信、生物、地震等领域。在实际情况中，很多数据都是非线性和非平稳的，而传统的信号处理方法，例如傅里叶变换、小波变换等，都是假设数据信号是线性和平稳的。1998年由NASA的Norden E Huang等人提出了一个广泛的信号分析技术，称为Hilbert-Huang变换（简称HHT）。它通过结合经验模态分解（EMD）和Hilbert频谱分析来适应非线性与非平稳数据的处理。
　　R提供了计算、绘图的计算机语言和环境，它提供了丰富的统计功能和绘图技术。更为重要的是R作为开源软件，围绕它已经建立起了开源生态圈，不同领域专家（包括数学、金融、通信、生物等）都利用R进行相关工作，并将其工作成果上传至开源社区共享。如今对于Hilbert-Huang变换的实现大多基于MATLAB，由此基于R实现Hilbert-Huang变换为R的使用者提供了有力的工具。
　　本文首先介绍Hilbert-Huang变换的基本原理，并分析实现Hilbert-Huang变换的过程和结构；其次说明如何使用R语言相关工具包实现Hilbert-Huang变换；最后，使用金融高频数据验证实现代码并给出相应结果。
　　2 Hilbert-Huang变换（HHT）
　　Hilbert-Huang变换是一种时间频率技术，它包括两个部分。第一部分，通过经验模态分解（EMD）把信号分解为固有模态函数（IMF），提取出信号中的尺度特征。第二部分，针对每个IMF进行Hilbert变换产生时间频率表示。
　　固有模态函数（IMF）是一个满足以下两个条件的函数：
　　1）在整个数据集中，极值数量和过零点数据必须相等或相差1；
　　2）在任何一点，由极大值和极小值定义的包络线的均值等于0。
　　大部分数据序列都不满足以上标准。在任意时间点，数据可能包含多个振荡模式。这就使简单的Hilbert变换不能提供给定数据频率的完整描述。因此，首先必须通过EMD程序从数据中分解振荡模式，获得固有模态函数（IMF），基本算法如下：
　　1）定义信号[x（t）]中所有的极值；
　　2）通过对最小值和最大值进行插值（一般通过三次样条插值），获得下包络线[xL（t）]和上包络线[xU（t）]；
　　3）计算上下包络线的均值[e（t）=xL（t）+xU（t）2]；
　　4）分解内在振荡模式[d（t）=x（t）-e（t）]；
　　5）让[d（t）]作为一个新信号[x（t）]。重复上述步骤，直到[d（t）]变为符合停止准则的零均值过程。
　　6）一旦得到了零均值[d（t）]，它就可以被指定为第一个固有模态函数IMF 1，c1。
　　7）IMF 1从原始信号中减去，剩余部分被用作新信号[x（t）]。重复筛选过程得到IMF 2，c2。
　　8）重复上述步骤，得到c3，c4等等。这个过程直到剩余部分是一个只有最小值和最大值的单调函数。
　　获取到IMF分量后，就可以通过对每个分量应用Hilbert变换（HT）来获取瞬时频率。对于信号[x（t）]的HT被定义为[H[x（t）]=x×1πt=y]，或使用卷积定义[y（t）=1πP-∞+∞x（τ）t-τdτ]，
　　其中[P]为柯西主值。通过[y（t）]定义分析信号[z（t）=x（t）+iy（t）]或[z（t）=a（t）eiθ（t）]，其中[a（t）=x2（t）+y2（t）]，[θ（t）=arctany（t）x（t）]。
　　瞬时频率可以通过瞬时相位变化来定义[ω=dθ（t）d（t）]，这就是HHT的第二部分，即Hilbert频谱分析（HSA），它必须通过以下形式记录信号[x（t）=j=1naj（t）eiωj（t）dt]。
　　由此，获取到作为时间函数的振幅和瞬时频率，这个时间频率的表示就是Hilbert频谱（HS）。为了获取每个频率值对总振幅（能量）的贡献，可以计算出边际谱，定义如下：[h（ω）=0TH（ω，t）dt]。
　　3 使用R实现HHT
　　R的开源性为广大使用者提供了大量可用的包，其中由Kim和Oh在2008年开发的EMD包提供了HHT的关键实现。现在我们讨论如何使用EMD包实现EMD和Hilbert频谱。
　　3.1 经验模态分解（EMD）
　　Hilbert-Huang变换的第一部分是通过经验模态分解（EMD）把信号分解为固有模态函数（IMF），这个过程中需要取得信号中的极值和过零点。EMD包提供了extrema函数来求取。extrema函数返回一个list，其中包含minindex、maxindex、nextreme、cross及ncross，分别代表局部最小值的时间索引矩阵、局部最大值的时间索引矩阵、极值数量、过零点的时间索引矩阵及过零点的数量。通过extremadata <- extrema（data）代码计算data序列的相关数据。
　　在前面提到了Hilbert-Huang变换的基本步骤，其中前6步获取固有模态函数（IMF）。EMD包中提供了extractimf函数，它从一个给定信号中分解所有IMF。extractimf函数提供了几个参数，包括residue、tt、tol、、stoprule、boundary、sm、spar、check，它们的含义分别为tt时间中观察到的数据或信号、观察的时间、停止筛选过程的容忍值、最大分解次数、停止筛选的条件、边界条件、包络线的构建方式、构建包络线参数及是否显示每次循环。其中大部分参数都有默认值，可以简单的通过imfData <- extractimf（data， time，check=TRUE）代码提取IMF。extractimf函数返回imf、r（分解后的剩余部分）及循环的次数。
　在extractimf函数中，有几个参数需要详细讨论一下。
　　首先，筛选过程的停止条件，stoprule参数可以选择"type1"或"type2"。"type1"指定当所有t时间满足扩展原有信号，添加部分通过靠近边界的极值形成的固有模式构建。Zeng and He （2004）提出两种扩展信号的方式，它们分别通过对称和自反方式添加信号，被称为偶扩展和奇扩展。可以通过指定boundary参数为" wave "、" symmetric "、" periodic "、" evenodd "或" none "来改变边界调整的方式。
　　通过整合extrema和extractimf函数就可以实现EMD，EMD包也提供了emd函数直接实现整个过程。emd函数的参数和extractimf类似，emd函数多了和两个参数，分别指定IMF的最大数量及是否显示每个IMF。
　　3.2 Hilbert频谱
　　实际工作中大部分信号都是非平稳的，也就是说信号的能量和振幅都不断在变化。得到IMF之后可以通过 Hilbert频谱来获取瞬时频率，以此来获取信号在各时间点能量和振幅的变化。EMD包中提供了hilbertspec函数计算信号的振幅和瞬时能量。hilbertspec函数有两个参数：xt和tt，分别指定信号矩阵和时间序列。hilbertspec函数返回一个振幅和瞬时频率的矩阵。EMD包还提供了spectrogram函数，它产生一张时间序列的振幅和瞬时频率图。该图的x轴表示时间，y轴表示瞬时频率，图片中每个点代表在某时间点某个频率的振幅。spectrogram函数的参数包括：amplitude，freq，tt，multi，nlevel，size，分别代表振幅（由hilbertspec函数获取）、频率（由hilbertspec函数获取）、时间序列、是否为多个信号、图片中点的颜色种类以及图片大小。
　　4 数据实验
　　本节采用真实数据进行Hilbert-Huang变换，我们使用2004年上证综指1月份5分钟高频交易数据为实验数据。图1中第一行显示了原始数据。
　　首先通过EMD函数分解所有IMF，图1中从第二行开始显示了6次筛选获取的IMF以及最后筛选完成后剩余的r。
　　其次，使用hilbertspec函数计算每个imf的振幅和瞬时能量。
　　5 结论
　　通过R的EMD开源包可以快速的实现Hilbert-Huang变换，为非统计专业领域的专家提供了使用Hilbert-Huang变换的有效工具。该文讨论了Hilbert-Huang变换的基本实现原理，并介绍了EMD包提供的基本函数的使用，最后通过真实的金融数据进行实验，并给出了实验结果。
　　参考文献：
　　[1] Norden E. Huang & Samuel S. P. Shen， The Hilbert-Huang Transform and Its Applications， World Scientific Publishing Co. Pte. Ltd，2005.
　　[2] N. E. Huang， Z. Shen， S. R. Long， M. C. Wu， H. ， Q. Zheng， N.-C. Yen， C. C. Tung， and H.H. Liu. The Empirical Mode Decomposition and Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis. Proceedings of the Royal Society London A.， 1998（454）：903-995.
　　[3] N. E. Huang， M. C. Wu， S. R. Long， S. Shen， W. Qu， P. Gloerson， and K. L. Fan. A Confidence Limit for the Empirical Mode Decomposition and Hilbert Spectral Analysis. Proceedings of the Royal Society London A.， 2003（31）：417-457.