CN105760347B - 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 - Google Patents
一种基于数据/极值联合对称延拓的hht端点效应抑制方法 Download PDFInfo
- Publication number
- CN105760347B CN105760347B CN201610077796.XA CN201610077796A CN105760347B CN 105760347 B CN105760347 B CN 105760347B CN 201610077796 A CN201610077796 A CN 201610077796A CN 105760347 B CN105760347 B CN 105760347B
- Authority
- CN
- China
- Prior art keywords
- symmetric
- signal
- point
- data
- continuation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Complex Calculations (AREA)
Abstract
本发明提供一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，包括如下步骤：步骤1、对原始信号进行端部数据对称延拓；步骤2、极值点对称延拓；步骤3、获取联合延拓后信号的包络线；步骤4、EMD分解，截取原始信号时域范围内的IMF分量；步骤5、将包含延拓部分的IMF分量进行Hilbert变换，截取原始信号时域范围内的瞬时频率，解决EMD分解中固有的端点飞翼问题以及Hilbert变换过程中的端点发散问题。
Description
技术领域
本发明涉及一种基于数据/极值联合对称延拓的HHT端点效应抑制方法。
背景技术
传统的信号分析与处理大多数建立在傅立叶变换(Fourier Transform，简称FT)的基础上。但是，傅立叶变换只能处理全域波，并要求被分析系统必须是线性的、信号必须是平稳或者是严格周期的，这就大大限制了傅立叶变换在非线性系统、非平稳信号下的应用。并且，尽管傅立叶变换在频域或时域内能获得较高的分辨率，但它不能在这两个领域内同时满足，无法表征出信号中的某一频率分量出现的时间点，以及分量随着时间的变换趋势，因而无法满足非平稳信号分析处理的需要。后来出现的短时傅立叶变换、Wigner-Ville分布、小波变换，虽然在一定程度上能处理非平稳信号，但严格说来，它们仍属于全局分析方法的范畴，对信号的局部没有自适应性，不能精确表述频率随时间的变化，即不能做时频局部分析。
为了进一步分析处理非线性、非平稳信号，1998年，Norden E.Huang提出了一种新的信号处理方法：经验模态分解法(Empirical Mode Decomposition，简称EMD)，并引入了基于Hilbert变换的Hilbert谱的概念和Hilbert谱分析的方法，美国宇航局(NASA)将该方法命名为Hilbert-Huang Transform，简称HHT。EMD的本质是对非平稳信号进行平稳化处理，其结果是将非平稳信号中不同尺度的波动或趋势逐级分解开来，产生一系列具有不同特征尺度的数据系列，每个序列称为一个固有模态函数(Intrinsic Mode Function，简称IMF)。固有模态函数应满足两个条件：(1)信号的极值点的数目与过零点的数目相等或至多相差一个；(2)信号的局部极大值和局部极小值定义的包络均值为零。因此基于这些分量进行Hilbert变换得到的结果能够反映真实的物理过程，即信号能量在空间(或时间)各种尺度上的分布规律。
HHT方法处理信号的基本过程分为两部分，即经验模态分解(EMD)和Hilbert谱分析(Hilbert Spectral Analysis，简称HAS)；首先利用EMD方法将给定的信号分解为若干个固有模态函数(IMF)；然后将HAS作用在每一个IMF分量上，得到相应Hilbert谱；最后，汇总所有IMF分量的Hilbert谱获得原始非平稳信号的Hilbert谱。每一个有效的IMF分量应满足两个条件：(1)信号的极值点的数目与过零点的数目相等或至多相差一个；(2)信号的局部极大值和局部极小值定义的包络均值为零。
HHT在处理非平稳信号时存在以下两类端点效应问题：
1、EMD分解的端点飞翼现象
EMD分解中固有存在的端点效应直接影响到分解的有效性。在每一次的EMD筛分过程中，要根据信号x(t)的上、下极值点三次样条插值包络线来计算信号的局部平均值。由于信号是有限长的，信号左右两端点不能确定是局部极值点。如果直接在端点处弃值，那么进行三次样条插值时必然使信号的上下包络线在信号两端附近严重扭曲失真(如图1)，产生端点飞翼。特别是对于低频分量，由于其时间尺度大，极值间的距离大，这种端部边缘效应传播到信号的内部，严重影响了EMD分解的质量，使得分解出来的IMF分量没有实际的物理意义。
2、Hilbert变换的端点发散现象
对各个IMF分量进行Hilbert变换时，由于Hilbert变换是基于傅立叶变换实现的，而傅立叶变换会产生频谱泄露，出现严重的端点发散现象，这样得到的Hilbert谱并不会真实地反映原始信号的特征。(如图2)
针对上述EMD分解的端点飞翼现象，目前采用最多的是通过单次端点延拓的方法来加以抑制。主要方法分为两类：一类是延拓极值点，如极值点对称延拓、平行极值点延拓、多项式拟合延拓方法、极值平移法；另一类是延拓原信号数据，如采用神经网络延拓法、镜像闭合延拓法、AR模型延拓法、支持向量回归机法、边界波形匹配预测法、时变参数ARMA模型、最大相关波形延拓等。
针对Hilbert变换端点发散现象，目前主要通过神经网络、ARMA模型和支持向量回归机等数据延拓法先对IMF分量进行端点延拓，再进行Hilbert变换，并通过抛弃延拓的两端数据，使端点效应释放到原始数据的外端，保证原始数据范围内Hilbert变换的有效性。
现有技术缺点如下：
极值点延拓的方法由于不延拓信号本身数据点，仅以端点附近极值点的特性决定延拓极值，这对于宽频带、非平稳的复杂信号而言，显然不能准确反应信号的真实趋势。并且，由该方法所得各IMF分量之间正交性较差，Hilbert变换得到的瞬时频率端点发散严重。
神经网络预测算法需要训练的样本足够多，计算速度慢，耗时多，不易于实时信号的处理，而且具有局部极小点、过学习以及结构和类型的选择过分依赖于经验等固有的缺陷。
支持向量回归机理论复杂且不完善，比如，如何根据具体问题选择适当的内积函数，选择不同的参数(如损失函数，核函数，精度参数和惩罚参数等)及其不同的参数对延拓效果的影响等还有待进一步研究。
并且，现有方法多数通过2个阶段分别进行端点延拓实现EMD分解与Hilbert变换的端点效应抑制。
术语解释：
1、经验模态分解(Empirical Mode Decomposition，简称EMD)：将信号中不同尺度的波动或趋势逐级分解开来，产生一系列具有不同特征尺度的数据系列。
2、固有模态函数(Intrinsic Mode Function，简称IMF)：由EMD分解产生一系列具有不同特征尺度的数据系列，每个序列称为一个固有模态函数。固有模态函数满足两个条件：(1)信号的极值点的数目与过零点的数目相等或至多相差一个；(2)信号的局部极大值和局部极小值定义的包络均值为零。
3、Hilbert-Huang变换(Hilbert-Huang Transform，简称HHT)：指对非平稳信号进行EMD分解，进而对分解出的固有模态函数进行Hilbert变换，获得具有明确瞬时频率意义的Hilbert谱。
4、端部数据对称延拓：比较信号的端点值与端部极值，寻找合理对称点，向信号外侧镜像延长端部数据。
5、极值点对称延拓：比较信号的端点值与端部极值，寻找合理对称点，向信号外侧镜像添加端部极值。
6、数据/极值联合对称延拓：在端部数据对称延拓的基础上进一步进行极值点对称延拓。
发明内容
本发明要解决的技术问题，在于提供一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，解决EMD分解中固有的端点飞翼问题以及Hilbert变换过程中的端点发散问题。
本发明是这样实现的：一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，包括如下步骤：
步骤1、将原始信号进行端部数据对称延拓；
步骤2、极值点对称延拓；
步骤3、获取联合对称延拓后信号的包络线；
步骤4、EMD分解，截取原始信号时域范围内的IMF分量；
步骤5、将包含延拓部分的IMF分量进行Hilbert变换，截取原始信号时域内的瞬时频率。
进一步地，所述步骤1进一步具体为：设原始信号为x(t)，对x(t)执行如下操作：
步骤11、求x(t)的极大值点为U(1)，U(2)，…，U(n)，和x(t)的极小值点为L(1)，L(2)，…，L(m)；其中，m＝n，或丨m-n丨＝1，且m和n为整数；左端点记为x(0)，右端点记为x(end)；
步骤12、对原始信号的左端数据进行对称延拓：
a、当x(0)≥U(1)时，以x(0)为对称点，将x(0)～U(2)之间的信号数据向左对称延拓；
b、当x(0)≤L(1)时，以x(0)为对称点，将x(0)～L(2)之间的信号数据向左对称延拓；
c、当L(1)<x(0)<U(1)，且x(0)与U(1)相邻时，以U(1)为对称点，将U(1)～U(3)之间的信号数据向左对称延拓；
d、当L(1)<x(0)<U(1)，且x(0)与L(1)相邻时，以L(1)为对称点，将L(1)～L(3)之间的信号数据向左对称延拓；
步骤13、对原始信号的右端数据按步骤12同样的方法向右对称延拓；
步骤14、延拓后信号记为x1(t)；
所述步骤12以及步骤13不分先后顺序。
进一步地，步骤2进一步具体为：对经过端部数据对称延拓的信号x1(t)执行如下操作，将时间零点定于x1(t)左端点；
步骤21、求信号x1(t)的极大值点为U1(1)，U1(2)，…，U1(n)，极小值点为L1(1)，L1(2)，…，L1(m)；其中，m＝n，或者丨m-n丨＝1，且m和n为整数；其中左端点记为x1(0)，右端点记为x1(end)；
步骤22、对信号x1(t)的左端进行极值点对称延拓
a、当x1(0)≥U1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；
b、当x1(0)≤L1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；
c、当L1(1)<x1(0)<U1(1)，且x1(0)与U1(1)相邻时，以U1(1)为对称点，对称延拓极大值点U1(2)、U1(3)和极小值点L1(1)、L1(2)；
d、当L1(1)<x1(0)<U1(1)，且x1(0)与L1(1)相邻时，以L1(1)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(2)、L1(3)；
步骤23、对信号x1(t)的右端按步骤22同样的方法将极值点向右对称延拓。
进一步地，所述步骤3进一步具体为：
步骤31、确定信号x1(t)所有局部极值点(包括步骤22和23延拓的极值点)，用三次样条曲线连接所有的局部极大值和局部极小值，分别形成上包络线x1max(t)和下包络线x1min(t)。
进一步地，步骤4进一步具体为：
步骤41、求出上包络线x1max(t)与下包络线x1min(t)的平均值，记为m11(t)，将信号x1(t)去掉该平均值后得到新数据序列h11(t)：
m11(t)＝[x1max(t)+x1min(t)]/2
h11(t)＝x1(t)-m11(t)
步骤42、判断h11(t)是否满足IMF的两个条件，若不满足，则将h11(t)重复步骤21至41，直到新序列：
h1k(t)＝h1(k-1)(t)-m1k(t)
满足IMF的两个条件，得到第一个IMF分量c1(t)：
c1(t)＝h1k(t)
步骤43、从x1(t)中分离出c1(t)，得到剩余序列r1(t)：
r1(t)＝x1(t)-c1(t)
步骤44、将r1(t)作为新的序列，按照以上步骤21到步骤43，依次提取第n个IMF分量cn(t)；当残量rn(t)成为一个单调函数或小于某一预定值时，分解结束；
步骤45、截取步骤44中IMF分量ci(t)在原始信号时域范围内的数值，即为原始信号的IMF分量。
进一步地，所述步骤5进一步具体为：
将包含延拓部分的IMF分量ci(t)进行Hilbert变换得到
构造解析信号
得到相位函数
进一步取得瞬时频率截取原始信号时域范围内的瞬时频率。
本发明具有如下优点：
本发明无需在EMD分解和Hilbert变换两个阶段分别进行信号延拓以分别解决两阶段的端点效应问题，而是通过端部数据对称延拓，同时把两个阶段产生的端点效应释放到数据外端，一举解决两阶段的端点飞翼问题。
本发明在端点数据对称延拓的基础上，通过极值点对称延拓进一步控制三次样条包络曲线的走势，使原始数据的包络曲线更加准确。与现有的极值点对称延拓方法以及仅进行端部数据对称延拓的方法相比，本发明不仅具有更加良好的端点抑制效果，而且有利于改善信号的整体正交性。
与神经网络、ARMA模型和支持向量回归机的数据延拓等方法相比，本发明理论与算法简单，容易实现且计算速度快。
附图说明
下面参照附图结合实施例对本发明作进一步的说明。
图1为端点飞翼示意图。
图2为Hilbert变换端点发散示意图。
图3为信号x(t)左端数据对称延拓。
图4为信号x1(t)左端极值点对称延拓。
图5为获取联合延拓后信号的包络线。
图6为本发明一种基于数据/极值联合对称延拓的HHT端点效应抑制方法的详细流程图。
图7为简谐叠加信号。
图8为采用端点值作为极值点得到的IMF分量。
图9为采用极值点对称延拓得到的IMF分量。
图10为采用端部数据对称延拓得到的IMF分量。
图11为采用本发明方法得到的IMF分量。
图12为采用端点值作为极值点得到的瞬时频率。
图13为采用极值点对称延拓得到的瞬时频率。
图14为采用端部数据对称延拓得到的瞬时频率。
图15为采用本发明方法得到的瞬时频率。
具体实施方式
如图6所示，一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，包括如下步骤：
步骤1、对原始信号进行端部数据对称延拓：设原始信号为x(t)，对x(t)执行如下操作：
步骤11、求x(t)的极大值点为U(1)，U(2)，…，U(n)，和x(t)的极小值点为L(1)，L(2)，…，L(m)。其中，m＝n，或丨m-n丨＝1，且m和n为整数。左端点记为x(0)，右端点记为x(end)。
步骤12、对原始信号的左端数据进行对称延拓：
a、当x(0)≥U(1)时，以x(0)为对称点，将x(0)～U(2)之间的信号数据向左对称延拓。
b、当x(0)≤L(1)时，以x(0)为对称点，将x(0)～L(2)之间的信号数据向左对称延拓。
c、当L(1)<x(0)<U(1)，且x(0)与U(1)相邻时，以U(1)为对称点，将U(1)～U(3)之间的信号数据向左对称延拓。
d、当L(1)<x(0)<U(1)，且x(0)与L(1)相邻时，以L(1)为对称点，将L(1)～L(3)之间的信号数据向左对称延拓。
步骤13、对原始信号的右端数据按步骤12同样的方法向右对称延拓；
步骤14、延拓后信号记为x1(t)；所述步骤12以及步骤13不分先后顺序；
步骤2、极值点对称延拓：对经过端部数据对称延拓的信号x1(t)执行如下操作，将时间零点定于x1(t)左端点；
步骤21、求信号x1(t)的极大值点为U1(1)，U1(2)，…，U1(n)，极小值点为L1(1)，L1(2)，…，L1(m)；其中，m＝n，或者丨m-n丨＝1，且m和n为整数；其中左端点记为x1(0)，右端点记为x1(end)；
步骤22、对信号x1(t)的左端进行极值点对称延拓
a、当x1(0)≥U1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；
b、当x1(0)≤L1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；
c、当L1(1)<x1(0)<U1(1)，且x1(0)与U1(1)相邻时，以U1(1)为对称点，对称延拓极大值点U1(2)、U1(3)和极小值点L1(1)、L1(2)；
d、当L1(1)<x1(0)<U1(1)，且x1(0)与L1(1)相邻时，以L1(1)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(2)、L1(3)；
步骤23、对信号x1(t)的右端按步骤22同样的方法将极值点向右对称延拓；
步骤3、获取联合延拓后信号的包络线：
步骤31、确定信号x1(t)所有局部极值点(包括步骤22和23延拓的极值点)，用三次样条曲线连接所有的局部极大值和局部极小值，分别形成上包络线x1max(t)和下包络线x1min(t)；
步骤4、EMD分解，截取原始信号时域范围内的IMF分量：
步骤41、求出上包络线x1max(t)与下包络线x1min(t)的平均值，记为m11(t)，将信号x1(t)去掉该平均值后得到新数据序列h11(t)：
m11(t)＝[x1max(t)+x1min(t)]/2
h11(t)＝x1(t)-m11(t)
步骤42、判断h11(t)是否满足IMF的两个条件，若不满足，则将h11(t)重复步骤21和41，直到新序列：
h1k(t)＝h1(k-1)(t)-m1k(t)
满足IMF的两个条件，得到第一个IMF分量c1(t)：
c1(t)＝h1k(t)
步骤43、从x1(t)中分离出c1(t)，得到剩余序列r1(t)：
r1(t)＝x(t)-c1(t)
步骤44、将r1(t)作为一个新的序列，按照以上步骤21到步骤43，，依次提取第n个IMF分量cn(t)；当残量rn(t)成为一个单调函数或小于某一预定值时，分解结束。
步骤45、截取步骤44中IMF分量ci(t)在原始信号时域范围内的数值，即为原始信号的IMF分量。
步骤5、将包含延拓部分的IMF分量进行Hilbert变换，截取原始信号时域范围内的瞬时频率：
将包含延拓部分的IMF分量ci(t)进行Hilbert变换得到
构造解析信号
得到相位函数
进一步取得瞬时频率截取原始信号时域范围内的瞬时频率。
本发明可用于地震工程、地球物理探测、结构损伤侦探、机械工程、信息科学、海洋和大气科学、经济学、生态学、医学等领域所有涉及数据处理的科研和工程应用。
本发明一种具体实施方式如下：
4.1端部数据对称延拓
设原始信号为x(t)，对x(t)执行如下操作：
(1)求x(t)的极大值点为U(1)，U(2)，…，U(n)，和x(t)的极小值点为L(1)，L(2)，…，L(m)。其中，m＝n，或丨m-n丨＝1，且m和n为整数。左端点记为x(0)，右端点记为x(end)。
(2)对原始信号的左端数据进行对称延拓：
a、当x(0)≥U(1)时，以x(0)为对称点，将x(0)～U(2)之间的信号数据向左对称延拓。(如所图3中a示)
b、当x(0)≤L(1)时，以x(0)为对称点，将x(0)～L(2)之间的信号数据向左对称延拓。(如所图3中b示)；
c、当L(1)<x(0)<U(1)，且x(0)与U(1)相邻时，以U(1)为对称点，将U(1)～U(3)之间的信号数据向左对称延拓。(如所图3中c示)；
d、当L(1)<x(0)<U(1)，且x(0)与L(1)相邻时，以L(1)为对称点，将L(1)～L(3)之间的信号数据向左对称延拓。(如所图3中d示)；
(3)信号右端按同样方法向右对称延拓端部数据，延拓后信号记为x1(t)。
4.2极值点对称延拓
对经过端部数据对称延拓的信号x1(t)执行如下操作，为下文表述方便，将时间零点定于x1(t)左端点。
(1)求信号x1(t)的极大值点为U1(1)，U1(2)，…，U1(n)，极小值点为L1(1)，L1(2)，…，L1(m)。其中，m＝n，或丨m-n丨＝1。左端点记为x1(0)，右端点记为x1(end)。
(2)对信号x1(t)的左端进行极值点对称延拓
a、当x1(0)≥U1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；(如所图4中a示)
b、当x1(0)≤L1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；(如所图4中b示)
c、当L1(1)<x1(0)<U1(1)，且x1(0)与U1(1)相邻时，以U1(1)为对称点，对称延拓极大值点U1(2)、U1(3)和极小值点L1(1)、L1(2)；(如所图4中c示)
d、当L1(1)<x1(0)<U1(1)，且x1(0)与L1(1)相邻时，以L1(1)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(2)、L1(3)；(如所图4中d示)
(3)信号x1(t)的右端按步骤(2)同样方法将极值点向右对称延拓
4.3获取联合对称延拓后信号的包络线确定信号x1(t)所有局部极值点(包括步骤4.2延拓的极值点)，用三次样条曲线连接所有的局部极大值和局部极小值，分别形成上包络线x1max(t)和下包络线x1min(t)。(如图5所示)
4.4EMD分解
(1)求出上包络线x1max(t)和下包络线x1min(t)的平均值，记为m11(t)，将信号x1(t)去掉该平均值后得到新数据序列h11(t)：
m11(t)＝[x1max(t)+x1min(t)]/2
h11(t)＝x1(t)-m11(t)
(2)判断h11(t)是否满足IMF的两个条件，如果不满足，则将重复上述步骤4.2至4.4(1)处理过程，直到新数据序列
h1k(t)＝h1(k-1)(t)-m1k(t)
满足IMF的两个条件，得到第一个IMF分量c1(t)：
c1(t)＝h1k(t)
(3)从x1(t)中分离出c1(t)，得到剩余序列r1(t)：
r1(t)＝x1(t)-c1(t)
(4)将r1(t)作为一个新的序列，按照以上4.2至4.4(3)步骤，依次提取第二、三……直至第n个IMF分量cn(t)。当残量rn(t)成为一个单调函数或小于某一预定值时，分解结束。
(5)截取步骤(4)中IMF分量ci(t)在原始信号时域范围内的数值，即为原始信号的IMF分量。
4.5将包含延拓部分的IMF分量进行Hilbert变换，截取原始信号时域内的瞬时频率；
将包含延拓部分的IMF分量ci(t)进行Hilbert变换得到
构造解析信号
得到相位函数
进一步取得瞬时频率截取原始信号时域范围内的瞬时频率。
计算过程如图6所示。
以图7一个简谐叠加信号为例，按本发明方法进行HHT分析。已知信号：x(t)＝cos(2πt/50)+0.6cos(2πt/25)+cos(2πt/100)+0.5cos(2πt/200)，t∈[1,300]。图8为端点值作为极值点避免端点飞翼，其中虚线表示真实分量。图中IMF1分量端部信号明显收缩，最后低频分量IMF3的内部受到严重污染。
图9为仅对该信号进行极值点对称延拓处理(不执行步骤4.1)获得的IMF分量。IMF1与IMF2分量在端部与真实分量存在偏差，随着频率减小，低频分量IMF3的端点效应逐渐传播到信号内部。
图10为仅对该信号进行端部数据对称延拓(不执行步骤4.2)获得的IMF分量，IMF1与IMF2分量与真实分量吻合较好，但随着频率减小，低频分量IMF3的端点效应逐渐传播到信号内部，与真实分量产生较大误差。
图11为利用本发明提出的数据/极值联合对称延拓方法得到的分解结果，与图10相比，不仅在高频分量IMF1、IMF2与真实分量吻合较好，低频分量IMF3也几乎接近真实信号。说明其端点效应抑制效果更加良好。
从整体正交性指标IO值判断端点效应处理效果，IO值越小，正交化越好。从表1可得，本发明提出的数据/极值联合对称延拓处理法得到的IO值最小，端点值作为极值点的方法得到IO值最大，说明本发明提出的数据/极值联合延拓法的端点处理效果最优，支持向量回归机延拓次之。
通过运算时间的比较可知，支持向量回归机延拓的方法用时最长，远大于其他方法，本发明提出的方法运行时间仅0.312秒，具有很高的计算效率。
图12是上述叠加信号x(t)采用端点值作为极值点的方法进行HHT变换后的瞬时频率，由Hilbert变换产生的端点发散效应清晰可见。
图13是采用极值点对称延拓后进行HHT变换得到的瞬时频率，三个分量的边缘效应虽有所改善，但仍然没有完全得以抑制。
图14是采用端部数据对称延拓方法后进行HHT变换得到的瞬时频率，端点发散虽有显著的抑制，但仍在瞬时频率为0.005HZ的曲线右端出现相对幅度为45％的发散。
图15是采用本发明提出的数据/极值联合对称延拓方法得到的瞬时频率图，从图中可清楚看出三个分量频率分别对应叠加信号的三个频率成份0.04Hz、0.02Hz、0.005Hz。相比图14，瞬时频率为0.005HZ的曲线右端更加平稳，端点效应抑制得更加彻底。
表1各类延拓方法评价指标
虽然以上描述了本发明的具体实施方式，但是熟悉本技术领域的技术人员应当理解，我们所描述的具体的实施例只是说明性的，而不是用于对本发明的范围的限定，熟悉本领域的技术人员在依照本发明的精神所作的等效的修饰以及变化，都应当涵盖在本发明的权利要求所保护的范围内。
Claims (5)
1.一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，其特征在于：包括如下步骤：
步骤1、设原始信号为x(t)，对x(t)执行如下操作：
步骤11、求x(t)的极大值点为U(1)，U(2)，…，U(n)，和x(t)的极小值点为L(1)，L(2)，…，L(m)；其中，m＝n，或丨m-n丨＝1，且m和n为整数；左端点记为x(0)，右端点记为x(end)；
步骤12、对原始信号的左端数据进行对称延拓：
a、当x(0)≥U(1)时，以x(0)为对称点，将x(0)～U(2)之间的信号数据向左对称延拓；
b、当x(0)≤L(1)时，以x(0)为对称点，将x(0)～L(2)之间的信号数据向左对称延拓；
c、当L(1)<x(0)<U(1)，且x(0)与U(1)相邻时，以U(1)为对称点，将U(1)～U(3)之间的信号数据向左对称延拓；
d、当L(1)<x(0)<U(1)，且x(0)与L(1)相邻时，以L(1)为对称点，将L(1)～L(3)之间的信号数据向左对称延拓；
步骤13、对原始信号的右端数据按步骤12同样的方法向右对称延拓；
步骤14、延拓后信号记为x1(t)；
所述步骤12以及步骤13不分先后顺序；
步骤2、极值点对称延拓；
步骤3、获取联合对称延拓后信号的包络线；
步骤4、EMD分解，截取原始信号时域范围内的IMF分量；
步骤5、将包含延拓部分的IMF分量进行Hilbert变换，截取原始信号时域内的瞬时频率。
2.根据权利要求1所述的一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，其特征在于：步骤2进一步具体为：对经过端部数据对称延拓的信号x1(t)执行如下操作，将时间零点定于x1(t)左端点；
步骤21、求信号x1(t)的极大值点为U1(1)，U1(2)，…，U1(n)，极小值点为L1(1)，L1(2)，…，L1(m)；其中，m＝n，或者丨m-n丨＝1，且m和n为整数；其中左端点记为x1(0)，右端点记为x1(end)；
步骤22、对信号x1(t)的左端进行极值点对称延拓
a、当x1(0)≥U1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；
b、当x1(0)≤L1(1)时，以x1(0)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(1)、L1(2)；
c、当L1(1)<x1(0)<U1(1)，且x1(0)与U1(1)相邻时，以U1(1)为对称点，对称延拓极大值点U1(2)、U1(3)和极小值点L1(1)、L1(2)；
d、当L1(1)<x1(0)<U1(1)，且x1(0)与L1(1)相邻时，以L1(1)为对称点，对称延拓极大值点U1(1)、U1(2)和极小值点L1(2)、L1(3)；
步骤23、对信号x1(t)的右端按步骤22同样的方法将极值点向右对称延拓。
3.根据权利要求2所述的一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，其特征在于：所述步骤3进一步具体为：
步骤31、确定信号x1(t)所有局部极值点，所述局部极值点包括步骤22和23延拓的极值点，用三次样条曲线连接所有的局部极大值和局部极小值，分别形成上包络线x1max(t)和下包络线x1min(t)。
4.根据权利要求3所述的一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，其特征在于：步骤4进一步具体为：
步骤41、求出上包络线x1max(t)与下包络线x1min(t)的平均值，记为m11(t)，将信号x1(t)去掉该平均值后得到新数据序列h11(t)：
m11(t)＝[x1max(t)+x1min(t)]/2
h11(t)＝x1(t)-m11(t)
步骤42、判断h11(t)是否满足IMF的两个条件，若不满足，则将h11(t)重复步骤21至41，直到新序列：
h1k(t)＝h1(k-1)(t)-m1k(t)
满足IMF的两个条件，得到第一个IMF分量c1(t)：
c1(t)＝h1k(t)
步骤43、从x1(t)中分离出c1(t)，得到剩余序列r1(t)：
r1(t)＝x1(t)-c1(t)
步骤44、将r1(t)作为新的序列，按照以上步骤21到步骤43，依次提取第n个IMF分量cn(t)；当残量rn(t)成为一个单调函数或小于某一预定值时，分解结束；
步骤45、截取步骤44中IMF分量ci(t)在原始信号时域范围内的数值，即为原始信号的IMF分量。
5.根据权利要求4所述的一种基于数据/极值联合对称延拓的HHT端点效应抑制方法，其特征在于：所述步骤5进一步具体为：
将包含延拓部分的IMF分量ci(t)进行Hilbert变换得到
构造解析信号
得到相位函数
进一步取得瞬时频率截取原始信号时域范围内的瞬时频率。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610077796.XA CN105760347B (zh) | 2016-02-04 | 2016-02-04 | 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610077796.XA CN105760347B (zh) | 2016-02-04 | 2016-02-04 | 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105760347A CN105760347A (zh) | 2016-07-13 |
| CN105760347B true CN105760347B (zh) | 2019-04-23 |
Family
ID=56329954
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610077796.XA Active CN105760347B (zh) | 2016-02-04 | 2016-02-04 | 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105760347B (zh) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108009122B (zh) * | 2017-11-06 | 2021-01-05 | 天津大学 | 一种改进的hht方法 |
| CN108153710B (zh) * | 2017-12-15 | 2020-07-10 | 清华大学 | 基于静态小波变换系数最小化的多成分信号延拓方法 |
| CN108845230A (zh) * | 2018-06-22 | 2018-11-20 | 国网陕西省电力公司电力科学研究院 | 一种次同步振荡随机时变模态辨识方法 |
| CN109613609B (zh) * | 2019-01-16 | 2019-09-24 | 国家深海基地管理中心 | 一种基于变长度处理的组合信号分解方法 |
| CN111695226B (zh) * | 2019-03-12 | 2023-05-05 | 天津工业大学 | 一种用于信号分解的aclcd方法 |
| CN110514921B (zh) * | 2019-07-22 | 2022-07-26 | 华南理工大学 | 一种电力电子变换器非平稳信号中非线性现象的识别方法 |
| CN113362853B (zh) * | 2020-03-03 | 2022-07-01 | 东北大学秦皇岛分校 | 一种基于lstm网络emd端点效应抑制方法 |
| CN111444613B (zh) * | 2020-03-26 | 2023-09-08 | 云南电网有限责任公司电力科学研究院 | 基于改进emd电力系统谐波分析方法、装置及存储介质 |
| CN111580654A (zh) * | 2020-05-07 | 2020-08-25 | 重庆邮电大学 | 一种基于emd的脑电信号的短时特征提取方法 |
| CN111476220A (zh) * | 2020-06-03 | 2020-07-31 | 中国南方电网有限责任公司超高压输电公司大理局 | 一种换流阀空冷器故障定位方法 |
| CN114280526B (zh) * | 2022-03-03 | 2022-05-31 | 武汉格蓝若智能技术有限公司 | 一种电子式互感器校验仪数字微差溯源系统及方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102169476A (zh) * | 2011-04-14 | 2011-08-31 | 哈尔滨工业大学 | 一种基于灰色理论的希尔伯特-黄变换端点效应抑制方法 |
| WO2014144769A1 (en) * | 2013-03-15 | 2014-09-18 | Atrium Medical Corporation | Fluid analyzer and associated methods |
-
2016
- 2016-02-04 CN CN201610077796.XA patent/CN105760347B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102169476A (zh) * | 2011-04-14 | 2011-08-31 | 哈尔滨工业大学 | 一种基于灰色理论的希尔伯特-黄变换端点效应抑制方法 |
| WO2014144769A1 (en) * | 2013-03-15 | 2014-09-18 | Atrium Medical Corporation | Fluid analyzer and associated methods |
Non-Patent Citations (3)
| Title |
|---|
| 《Hilbert-Huang变换端点问题处理方法的分析与研究》;肖韵;《西南交通大学硕士学位论文》;20120401;第19,26-28页 |
| 《改进的迭代希尔伯特变换及滚动轴承故障诊断应用》;秦毅 等;《振动与冲击》;20130615;第32卷(第11期);第83-88页 |
| 《数值方法进展:从Fourier变换到Hilbert-Huang变换》;吴琛 等;《福建工程学院学报》;20151225;第13卷(第6期);第511-519页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105760347A (zh) | 2016-07-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105760347B (zh) | 一种基于数据/极值联合对称延拓的hht端点效应抑制方法 | |
| Shi et al. | A novel fractional wavelet transform and its applications | |
| Huang | An adaptive data analysis method for nonlinear and nonstationary time series: the empirical mode decomposition and Hilbert spectral analysis | |
| CN103956756A (zh) | 一种电力系统低频振荡模态辨识方法 | |
| CN107086566A (zh) | 基于广域信息的lmd互联电力系统低频振荡分析方法 | |
| CN101858938A (zh) | 基于自适应滤波原理的瞬时频率测量方法 | |
| CN109359633B (zh) | 基于希尔伯特-黄变换和小波脊线的信号联合分类方法 | |
| CN108334872A (zh) | 基于改进hht变换的特征提取方法 | |
| CN110726875A (zh) | 一种新能源柔性直流并网暂态谐波检测方法及系统 | |
| Hawley et al. | Some properties of an empirical mode type signal decomposition algorithm | |
| CN108761202B (zh) | 极点对称模态分解和希尔伯特变换相结合的谐波检测方法 | |
| CN108680782B (zh) | 基于极值点对称模式分解的电压闪变参数检测方法 | |
| Zhidong et al. | A new method for processing end effect in empirical mode decomposition | |
| CN109460614B (zh) | 基于瞬时带宽的信号时间-频率分解方法 | |
| Zhu et al. | An adaptive STFT using energy concentration optimization | |
| CN105466710B (zh) | 基于频域相似度的局部均值分解端点效应改进方法 | |
| CN109282841B (zh) | 一种无线无源声表面波传感器的超分辨率测量方法 | |
| Molla et al. | Hilbert spectrum in time-frequency representation of audio signals considering disjoint orthogonality | |
| CN108169560A (zh) | 一种分段正弦拟合分解方法 | |
| CN103577877A (zh) | 一种基于时频分析和bp神经网络的船舶运动预报方法 | |
| Pao et al. | Smoothing empirical mode decomposition: A patch to improve the decomposed accuracy | |
| Han | The analysis of signal based on the S-transform | |
| CN105405444B (zh) | 一种在Odd-DFT域对含噪正弦信号进行参数估计的方法 | |
| Zhou et al. | Morphological Filter‐Assisted Ensemble Empirical Mode Decomposition | |
| Zhang et al. | Determining the length of sliding window by using frequency decomposition |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |