经验模态分解算法中端点问题的处理(1).doc

x=0 30 60 90 120 150 180 210 240 270 300 330 360;y=-0.0167 -1.0927 -1.8725 -2.3586 -2.3061 -1.9576 -0.9574 -0.0080 0.8896 1.3877 1.1139 0.8517 -0.0167;fun=(a,t) a(1)+a(2)*sind(t+a(3) %matlab7.0以上版本，否则用inline%fun=inline(a(1)+a(2)*sind(t+a(3),a,t)a0=-0.5 -1.9 -0.079;a=nlinfit(x,y,fun,a0)t=0:5:360;yf=fun(a,t);plot(x,y,o,t,yf)结果：fun = (a,t) a(1)+a(2)*sind(t+a(3)a =-0.5239 -1.8995 -14.2382经验模态分解算法中端点问题的处理摘要：经验模态分解(EMD)方法就是对非线性、非平稳信号运用时间区域序列的上下包络线的均值得到瞬时平衡位置，将被分析信号分解成一组相互独立的稳态和线性的固有模态函数(IMF)数集。经验模态分解(EMD)方法是基于原始信号本事出发，经过筛选先把频率高的IMF分量分离出来，然后在分离频率较低的IMF分量。其实质就是利用时间特征尺度来获取原始信号数据中的振荡模态，本文对经验模态分解算法中端点问题的处理进行研究。关键词：经验模态分解 算法 端点 函数 经验模态分解(EMD)方法被提出后在各个领域普遍的应用，其具有直观、简单、自适应、完备性和正交性以及调制特性等一系列良好的特点。 (1)自适应性 经验模态分解(EMD)方法的自适应性表现为自适应生成基函数。在整个筛选分解过程中是根据原始信号自己的时间特征尺度实现的，不需要事先设立任何基函数。这与傅立叶变换和小波变换有着根本性的不同。傅立叶变换和小波变换需要事先设定谐波基函数和小波基函数，他们是先验性的。可以说在理论上，经验模态分解(EMD)方法适用任何信号的分解，其在对非线性、非平稳信号的处理上的优越性是其他时频分析方法无法比拟的。 经验模态分解(EMD)方法的自适应性还表现为自适应滤波特性。经验模态分解方法是基于原始信号本身出发，经过筛选先把频率高的IMF分量分离出来，然后在分离频率较低的工MF分量。这些不同频率成分以及带宽都是随原始信号的变化而改变的。因此，EMD方法可被视为是一组具有自适应性能的带通滤波器，它的截止频率和带宽均随着原始信号的变化而自动改变的，随着信号分析的目的改变而自动变化的。 这些尺度范围和频率成分均随着原始信号的变化而自动改变的，这样原始信号的特征可在不同分辨率下被表示，实现自适应多分辨率。 (2)正交性与完备性 所谓信号分解方法的完备性，是指可以从被分解后的信号的各个分量还原出原始信号的性质。经验模态分解(EMD)方法的本身就已经证明了其是完备的。可以得到证明EMD方法的完备性。且从EMD整个分解过程和结果都说明EMD方法的完备性。 所谓信号的正交性指的是被分解后的信号的各个分量之间相互正交的性质。例如频率不同的两个正弦信号它们是相互正交的。在EMD方法中，根据IMF的概念，每一个IMF分量应该是在局部应该是相互正交的。 此处的正交性是在局部意义上而言的。对于有一些特殊数据，有可能会出现两个相邻的分量在不同的时间段内含有相同的频率成分，因此在全局意义上不正交。由于一般在实际进行经验模态分解时采用的都是截取数据的有限长度，这样即使对于不同频率的纯正弦波形叠加信号进行分解也会有严重泄漏。泄漏的程度一般与数据的长度以及分解的结果是密切相关的。黄通过大量的实验数据证明EMD的泄漏一般小于1%;对于极短的数据为5%,与正弦型傅立叶分解在同一数量级上。据此可以认为EMD分解得到的各个IMF分量近似正交的。 (3) IMF分量的调制特性由固有模态函数的概念可知，对于任意信号被分解为有限工MF分量，这些分量可以是幅值和频率调制的。任何频率随时间的变化都可定义为频率调制。频率调制有两种概念:一是波间调制;二是波内调制。EMD分解得到的各个IMF分量不仅含有原始信号的非线性以及非平稳特性，而且工MF分量有波内调制特性，能用一个IMT表示由不同傅立叶频率描述的同一分量。1 经验模态分解和主成分分析1.1 经验模态分解经验模态分解算法的主要目的是将待分析信号分解为一系列表征时间尺度的IMF分量,要求IMF分量必须满足两个条件：IMF的极值点个数与过零点个数不超过1；由极大值点和极小值点确定的包络线均值为零.对信号进行EMD分解的步骤如下6：Step1 确定的所有极大值和极小值，分别对极大值点和极小值点进行三次样条插值,构造的上下包络线和，计算上下包络线的均值 ；Step2 计算和之间的差值, ； Step3 判断是否是一个IMF, 3.1) 如果符合IMF的定义条件，是一个IMF，则抽取作为第一个IMF分量,令，并求原信号与之间的差值， 。 3.2) 如果不是一个IMF,则将视为一个新的信号序列，重复步骤1和步骤2，求其包络均值及与间的差值。对重复上述过程次，直到符合IMF的定义条件，则令为的第一个IMF分量，并求原信号与之间的差值， ；Step4 将作为一个新的“原始”信号，重复步骤（step1-step3），抽取第2个内蕴模态函数分量，令；将作为一个新的“原始”信号，抽取第3个内蕴模态函数分量， ；以此类推，直到第次的余项满足终止条件，则停止迭代，的EMD分解完成。EMD分解结束后,原始信号可被表示为各IMFs和一个余项之和： 其中表示第个IMF分量.如果被零均值高斯白端点污染,则中所含端点仍近似服从零均值正态分布14,即可设 (1)其中表示没被污染的原始信号,表示所含端点,且.对于含随机端点的信号,先分解出的IMF分量通常对应于信号的高频端点,若去除几个先分解出的IMF,把剩余IMF进行重构,则可减弱信号的端点. 1.2 主成分分析主分量分解(PCA)是统计分析中常用的多通道数据分析方法15,被广泛应用于数据的降维和端点处理中.设道原始数据构成一个的数据阵,令,其中为的期望,记,则的协方差矩阵为 ,通过奇异值分解(SVD),协方差矩阵可被写为这里为对角矩阵,为的特征值,为特征值所对应的特征向量组成的正交矩阵.令 (2)称的各行为的主分量,它们在中依贡献率大小排序.对应的特征值与特征值总和的比称为该主分量的贡献率,表征该主分量代表原始信号能量的百分比.如果仅取前个主分量重构原始数据,则重构的数据为 端点信号在经PCA处理时,由于对应的主分量包含了信号中的大部分端点,在重构时直接丢掉；而代表了信号主要特征的主分量被保留,因而PCA可以有效的端点处理.2 利用3准则提取中的细节信息在最初的EMD端点算法中,通常认为第一层端点IMF全部由端点构成,并由此出发推导其他层IMF中所含端点的能量.但随着研究的深入,逐渐发现中仍含有一定量的信号细节信息12,13.对进行适当的处理,提取其所含的信号细节信息并加以保留,会提高端点效果,利用处理后的估计其余IMF中所含端点的能量也更准确.但由于先验知识很少,所以对进行处理是个难题.由EMD的研究可知,在中端点占绝大部分,而仅含有少量的信号细节信息,而且所含端点仍近似服从零均值正态分布,所以非常适合采用“3法则”进行细节信息提取16.由公式(1)可知, 满足加性端点模型 且.根据“3”法则,端点的分布满足即端点落在之间的概率为0.9973,而落在3之外的概率仅约为0.003.因此如果的值没有落在之内,则可认为中必然含显著误差,也即有必然含有信号信息,需要予以保留.利用“3法则”对进行细节信息提取可表示为： 其中表示从提取出的信号细节；端点方差采用文献17中提出的方法进行估计,即,这里表示的高频子带小波系数. 3 各层IMF中所含端点能量的估计利用“3法则”对进行细节信息提取后,可求出中所含端点的能量18： (3) 假设()中所含端点的能量为,则是未知的,由于信号和端点混杂在一起,因此一般情况下并不能求出.但通过含噪信号经EMD分解后的端点能量模型,可对进行近似计算.被白端点污染的信号经EMD分解后,如果第一层内蕴模态函数中所含端点的能量为,则中所含的端点的能量可由下式求出18, , (4)其中,.因此，求出中所含端点的能量后, 即可通过公式(4)估计中所含端点的能量.4 根据端点能量利用PCA去除()中的端点PCA是一种自适应的分解方法,信号经PCA分解后各主分量间互不相关,而且按照贡献率选择合适主分量重构后,能有效端点处理并保留信号绝大部分的主特征信息. ()经PCA分解后,信号和端点能够被有效分离,如果想较好去除中端点,必须要选择合适个数的主分量进行重构,通常根据前个主分量的累计贡献率来确定保留的主成分分量的个数.但累积贡献率的选择并不是一个简单的事情：1)累计贡献率取得太大,会残留较多的端点,导致端点不能完整的去除；累计贡献率取得太小,又会损失较多的信号细节信息；2) 每一层端点项IMF中所含端点的强度并不相同,因此在对不同层的IMF进行PCA端点时,累计贡献率不能取固定的值.为了更好地端点处理,本文根据中所含端点能量的比例,提出了一种自适应确定累计贡献率的方法.由公式(1)可知 ,().为了表示方便,设,为的协方差矩阵,其中表示的长度,显然所含端点与所含端点相同.假设经PCA分解后的主分量为,如果选择前 个主分量进行重构以去除中的端点,则可得到端点后的信号 (5)此时从中删除的端点为 (6 )设和的能量分别为、.在PCA端点中,通常认为前几项主分量包含了信号的主要特征信息,而比较靠后的主分量主要由端点构成,而且按照主分量对总能量的贡献率选择应保留的主分量个数.因此在利用PCA对端点时,如果选择合适的,使删除的端点的能量与本身所含的端点能量相同,也即使得 则可认为中的端点基本被全部去除,达到了较好的端点效果.上式等价于 (7)由公式(2)和的正交性可知： ,所以信号的能量为： (8)而所删除的端点的能量可表示为 (9)由公式（8）和（9）可知,被删除端点的能量与的能量之比为.所以,为了使中的端点被完整去除,应选择合适的使成立 .但在选择时,很难保证使得恰好成立,本文中对按照以下方法进行取值：如果存在使得式(10)成立, 则令, (10)应保留的主分量个数确定后,根据公式(5)可求出端点后的信号,因为,所以端点后的值为：.本文所提出的基于PCA的EMD端点算法具体步骤为：Step1 对信号进行EMD分解,设分解后的IMF为, ,余项为；Step2 对采用“3法则”提取信号细节信息,设提取的细节信息为； Step3 根据公式(3)求所含端点的能量,并利用公式(4)估计()中所含端点的能量；Step4 对()进行PCA分解,根据公式(10)选择合适个数的主分量进行重构端点,设端点后的值为；Step5 累加全部()和余项,得到端点后信号.5 实验结果与分析当一维信号长度较大时,其协方差矩阵的规模较大,直接进行PCA变换运算量很大.为了降低运算复杂度,本文按照文献19中所提出的嵌入方法首先将一维信号转换为多维信号后,再进行PCA分解,此时分解结果的特征和性质均保持不变.为了分析所提出算法的端点性能,分别对模拟信号和真实信号进行端点实验.模拟信号利用Matlab中的wnoise函数生成,分别生成“Blocks”,“Bumps”,“Heavy sine”和“Doppler”等四类具有典型特征的测试样本；真实信号选择来自Bell实验室的一段电力系统信号. 为了对比端点效果,对端点信号分别采用基于Bayesian阈值的小波端点法(Bayesian-Wavelet)、基于模态单元的EMD阈值法(Mode-EMD)和本文提出的基于PCA的EMD端点法 (PCA-EMD) 进行端点. 在利用Bayesian-Wavelet法进行端点时，小波基选用“db8”小波,分解层数取为10；在Bayesian-Wavelet和Mode-EMD的端点中,均采用硬阈值法. 本文采用均方误差MSE和信噪比SNR来评估算法的性能：信噪比越大,均方误差越小,表明端点效果越好.5.1 对模拟信号的端点在对模拟信号的实验中,首先利用wnoise生成信噪比分别为SNR=0,5dB,10,15,20dB的测试信号,信号长度取为4096. 图1是SNR=5dB时,测试样本“Doppler”及不同方法端点后的实验结果,在图1(c)-(f)中,虚线为原始信号,实线为端点后信号.对3种方法端点后的结果分别计算MSE和SNR(见表1).可以看出,在当SNR=5dB时,“Doppler”信号经PCA-EMD方法端点后效果相对较好,与Bayesian-Wavlet算法相比,SNR提高了2.509,MSE约减小了0.053；与Mode-EMD算法先比,SNR约提高了1.529,MSE约降低了0.017.不同信噪比的模拟信号经三种方法端点后的MSE和SNR如表1所示,通过比较可知,本文提出的PCA-EMD方法总体端点效果要优于Bayesian-Wavelet方法和Mode-EMD方法,端点后的信号更接近原始信号；但当信噪比增大时,PCA-EMD与Mode-EMD之间的差距在逐渐减小,当信噪比增加到20dB时,本文方法与Mode-EMD算法的端点结果已非常接近.(a) (b) (c) (d) (e) Fig.1 Doppler signal de-noising results comparison (a)Original Doppler signal (b)noisy Doppler signal (d) Bayesian Wavelet de-noising result (d)Mode-EMD de-noising result (e) PCA-EMD de-noising result图1 Doppler信号的端点结果比较 (a)原始Doppler信号 (b)含噪Doppler信号 (c) Bayesian Wavelet端点 (d) Mode-EMD端点 (e)PCA-EMD端点Table 1 Results of experiments using noisy simulated signals表1 模拟信号的端点结果BlocksMethodsSNR/Variance0 dB5 dB10 dB15 dB20 dBPCA-EMD14.529/0.213816.820/0.135220.761/0.054522.959/0.027423.905/0.0257Mode-EMD14.149/0.233316.571/0.133620.390/0.055422.833/0.031623.225/0.0289Bayesian Wavelet13.895/0.247416.232/0.144419.822/0.063222.232/0.032923.537/0.0285BumpsMethodsSNR/Variance0 dB5 dB10 dB15 dB20 dBPCA-EMD14.061/0.020417.044/0.010219.619/0.005719.794/0.005117.291/0.0097Mode-EMD13.366/0.022416.643/0.011319.316/0.006119.152/0.005817.220/0.0099Bayesian Wavelet12.001/0.032815.943/0.013218.857/0.006818.939/0.006317.197/0.0097Heavy SineHeavy SineMethodsSNR/Variance0 dB5 dB10 dB15 dB20 dBPCA-EMD19.735/0.101322.124/0.058427.544/0.016831.206/0.008534.194/0.0036Mode-EMD18.612/0.131121.946/0.060927.197/0.018230.337/0.008834.142/0.0037Bayesian Wavelet18.414/0.137321.358/0.069726.863/0.019630.479/0.009233.072/0.0047DopplerMethodsSNR/Variance0 dB5 dB10 dB15 dB20 dBPCA-EMD18.612/0.012421.549/0.006126.823/0.022429.918/0.013931.513/0.0024Mode-EMD16.426/0.020720.020/0.007825.419/0.025627.355/0.015830.932/0.0022Bayesian Wavelet15.819/0.022519.040/0.011424.684/0.029126.539/0.015028.212/0.00415.2 对实际电力信号的端点图2(a)是在有端点干扰的环境下采集到的一段电力系统信号,对该信号分别利用三种方法进行端点,端点后的结果如图2(b)-(d)所示.从图2可以看出,PCA-EMD端点后信号的细节和突变部分保持较好,与Bayesian-Wavelet和Mode-EMD端点结果相比,端点明显有所减少.分别计算各端点信号的MSE和SNR,结果如下：Bayesian-wavelet端点后,MSE=81.8073,SNR=32.023；Mode-EMD端点后,MSE=73.0587, SNR=32.5142,；PCA-EMD端点后,MSE=70.9869,SNR=32.6391.可以看出,PCA-EMD方法端点后的MSE最小,而SNR最大.这表明电力信号经PCA-EMD方法端点后,信号细节的保留程度和信号的复原程度都更好一些.可见,对实际信号进行端点时,PCA-EMD方法与Bayesian-Wavelet和Mode-EMD方法相比,端点能力也有一定程度的提高,是一种比较有效的端点方法.实验在matlab7.8.0环境下进行, 计算机内存为2G，cpu主频为3.06GHz， EMD采用flandrin 提供的pack_EMD (http:/perso.ens-lyon.fr/patrick.flandrin/emd.html)。统计Doppler，Heavy Sine， Bumps，Blocks以及实际电力信号的端点时间开销：本文方法的平均时间开销约为0.032943秒；Mode-EMD算法的平均时间开销约为0.029326秒。本文方法的时间花费比Mode-EMD算法约高12%，运算量没有明显增加。(a) (b) (c) (d) Fig.2 Electric signal de-noising results comparison (a) original electric signal (b)Bayesian Wavelet de-noising result (c)Mode-EMD de-noising result (d)PCA-EMD de-noising result 图2 电力信号的端点结果比较 (a)原始电力信号 (b)Bayesian Wavelet小波端点 (c)Mode-EMD端点 (d)PCA-EMD端点6 结束语为了提高EMD的端点效果,本文提出了一种基于PCA的EMD端点新方法.在所提出的方法中,首先利用“3法则”对第一层IMF进行信号细节信息提取,并计算各层IMF中所含端点的能量；然后利用PCA对各IMF进行端点,端点时根据每层IMF中所含端点的能量自适应的确定应保留的主分量的个数. 为了验证所提算法的性能,选择了基于Bayesian阈值的小波方法和基于模态单元的EMD阈值方法进行比较.首先采用含不同强度端点的模拟信号进行实验,实验结果表明,所提算法的端点效果整体优于小波方法和基于模态的EMD方法.而对含噪电力信号进行端点实例应用,也进一步验证了本文方法的有效性.参考文献:D. L. Donoho. De-noising by soft-thresholding. IEEE Transaction on Information Theory, 1995,41(3):613627.A. Pizurica and W. Philips. Estimating the probability of the presence of a signal of interest in multiresolution single and multiband image denoising. IEEE Transaction on Image Processing, 2006, 15(3):654665.J. Portilla, V. Strela, M. J. Wainwright, and E. P. Simoncelli. Image denoising using scale mixtures of gaussians in the wavelet domain. IEEE Transaction on Image Processing, 2003, 12(11):13381351.L. Zhang, P. Bao, and X. Wu. Multiscale LMMSE-based image denoising with optimal wavelet selection. IEEE Trans. Circuits Syst.Video Technol., 2005,15(4):469481.Zhang Hua, Chen Xiaohong, Yang Haiyan. Optimistic wavelet basis selection in seismic signal noise elimination. OGP, 2011, 46(1):70-75.Huang N E, Shen Z, Long S R, etc. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysisJ. Proc. of the Royal Society of London,1998, A454: 903995.Xu Xiaogang, Xu Guanlei, Wang Xiaotong, etc. Empirical mode decomposition and its application. Acta electronica Sinica, 2009, 37(3):581-585.Boudraa A, Cexus J. EMD-based signal filteringJ. IEEE Transaction on Instrumentation and Measurement, 2007, 56(6):2196-2202.Gan Yu, Sui Lifen. De-noising method for gyro signal based on EMD. Acta Geodaetica et Cartographica Sinica,2011,40(6):745-750.Olufemi A, Vladimir A, Auroop R. Empirical mode decomposition technique with conditional mutual information for denoising operational sensor dataJ. IEEE Sensors Journal, 2011, 11(10):2565-2575.Zhao Zhihong, Yang Shaopu, Shen Yongjun. Improved EMD based on de-noising method. Journal of Vibration and Shock, 2009,28(12):35-37.Qu Congshan, Lu Tingzhen, Tan Ying. A modified empirical mode decomposition