CN108009122B - 一种改进的hht方法 - Google Patents

一种改进的hht方法 Download PDF

Info

Publication number
CN108009122B
CN108009122B CN201711081152.9A CN201711081152A CN108009122B CN 108009122 B CN108009122 B CN 108009122B CN 201711081152 A CN201711081152 A CN 201711081152A CN 108009122 B CN108009122 B CN 108009122B
Authority
CN
China
Prior art keywords
imf
signal
imf component
component
value
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
Application number
CN201711081152.9A
Other languages
English (en)
Other versions
CN108009122A (zh
Inventor
张涛
丁碧云
潘希豪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN201711081152.9A priority Critical patent/CN108009122B/zh
Publication of CN108009122A publication Critical patent/CN108009122A/zh
Application granted granted Critical
Publication of CN108009122B publication Critical patent/CN108009122B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Compression, Expansion, Code Conversion, And Decoders (AREA)
  • Complex Calculations (AREA)

Abstract

一种改进的HHT方法：对原始信号挑选异常的极值点提取间歇信号，运用AR包络预测方法弥补异常的极值点，得到修复后的信号，运用Sliding‑fastBSpline‑EMD方法分解信号，得到所有IMF分量；根据得到的所有IMF分量与原始信号，计算每一个IMF分量的相关系数、互信息量和能量比值三组指标，并将每一个IMF分量的三组指标与三组指标相对应的设定阈值进行比较，根据比较结果对该IMF分量进行筛选，最后得到能够真实反映信号特征的一组IMF分量；对该组IMF分量作基于AR模型包络预测的归一化希尔伯特变换，根据希尔伯特变换得到IMF分量的瞬时相位，对瞬时相位求导，最终得到信号的瞬时频率。本发明具有计算量小、准确的特点能够得到更准确的信号特征。

Description

一种改进的HHT方法
技术领域
本发明涉及一种HHT方法。特别是涉及一种用于非线性非平稳信号时频分析的改进的HHT方法。
背景技术
1、HHT简介
希尔伯特-黄变换(Hilbert-Huang Transform，HHT)是N.E.Huang于1998年提出的一种专门针对非线性、非平稳信号处理的新型时频分析方法。它是一种适合分析非平稳过程的信号处理方法。该方法主要分成两部分，首先对信号进行经验模式分解((EmpiricalMode Decomposition，EMD)得到一系列的本征模式函数(Intrinsic Mode Function，IMF)，然后对本征模式函数进行Hilbert变换。
其中每个固有模态函数必须满足以下2个条件：
(1)在整个数据范围内，极值点与过零点的数目必须相等或者最多相差一个；
(2)在任意一点处，所有极大值点形成的上包络线和所有极小值点形成的下包络线的平均值始终为零。
其中EMD分解的主要步骤如下：
(1)对于信号x(t)，找出信号的所有局部极大值、极小值点；
(2)对这些极值点进行三次样条插值，得到由所有局部极大值点构成的上包络线和所有局部极小值点构成的下包络线，分别记为u(t)、v(t)；
(3)上下包络线的均值为
Figure GDA0002715474440000011
(4)令h(t)＝x(t)-m(t)，检验h(t)是否满足IMF的两个条件，如果满足，那么h(t)为分解出的第一个IMF；如果不满足，以h(t)为输入继续进行前面的步骤，直至得到第一个IMF，并将其记为c1(t)；
(5)将r1(t)＝x(t)-c1(t)作为新的待分析信号，重复步骤(1)～(4)得到c2(t)，此时记r2(t)＝r1(t)-c2(t)，重复上述步骤直到得到余项rm(t)是一个单调信号或者其值小于某个预先给定的阂值，分解结束。
最终我们得到m个IMF分量c1(t),c2(t),…,cm(t)，以及余项rm(t)，那么原始信号就可以表示为
Figure GDA0002715474440000012
尽管HHT方法在处理非线性、非稳态信号方面有很大的优势，但还处在发展阶段，在理论的严密性和方法的实用性方面还有很多工作要做。HHT目前存在模态混叠、端点效应、包络拟合误差和瞬时频率误差问题。
2、消除虚假IMF分量的问题
由于EMD分解过程中存在着包络拟合、模态混叠、端点效应和过分分解等缺点，在分解中会产生一些虚假的IMF分量，不足以反映信号本身的特征，在实际应用中，需要剔除这些虚假的IMF分量。真实IMF分量与原始信号的关联性比伪分量与原始信号的关联性大，且真实IMF分量占的比例较虚伪假IMF分量要大。目前剔除虚假IMF较常用的方法有相关系数法、互信息量法和能量比值法，表达式如下：
1)互相关系数：
Figure GDA0002715474440000021
其中ri(t)是第i个IMF分量ci(t)与原始信号x(t)的互相关系数，m是信号经过EMD分解得到的所有IMF分量的个数；
2)互信息量：
对于信号的第i个IMF分量ci(t)与原始信号x(t)
Figure GDA0002715474440000022
其中p(ci)和p(x)分别是第i个IMF分量ci(t)与原始信号x(t)的边缘概率分布，p(ci,x)是第i个IMF分量ci(t)与原始信号x(t)的联合概率分布。
3)能量比值：
Figure GDA0002715474440000023
其中，N表示一个信号中采样点的个数，Eni是第i个IMF分量ci(t)的总能量，En原始信号x(t)的总能量，si表示第i个IMF分量的总能量占原始信号总能量的比值，通过si的值的大小来判断第i个IMF分量所占能量是否是虚假分量。
3、瞬时频率的误差问题
由于Hilbert变换有一定的适用范围，对IMF分量进行Hilbert变换得到的瞬时频率结果中存在误差，若使用Hilbert变换求出有意义的瞬时频率，信号所需具备的条件受到Bedrosian定理和Nuttall定理的限制，这会导致所得的瞬时频率存在许多异常值。
针对这个问题，归一化Hilbert变换(NHT)可以明显降低常规Hilbert变换时出现的瞬时频率误差问题。常规的NHT变换的步骤如下：
1)对每一个待处理的IMF分量ci(t)求绝对值，记为|ci(t)|，找到|ci(t)|中所有的局部极大值。
2)使用三次样条插值方法对所有的局部极大值进行拟合，同样在进行拟合时，通过AR模型对的端点进行延拓，最后拟合得到经验包络信号e1(t)。
3)利用经验包络信号e1(t)对IMF分量ci(t)进行归一化处理：
y1(t)＝ci(t)/e1(t)
然后，对y1(t)进行检验，若y1(t)≤1，则结束归一化过程；若y1(t)不满足上面的条件，则继续重复前面的步骤，直到最后的结果满足y1(t)≤1的条件为止或循环次数为预设的终止次数m为止，一般取m＝4。
其中
Figure GDA0002715474440000031
此时，yj(t)即为IMF分量ci(t)的调频部分，而ci(t)的调幅部分ae(t)可用下式表示：
ae(t)＝e1(t)e2(t)…ej(t)
根据上述的归一化步骤，可以把IMF分量ci(t)用下式的形式表达：
ci(t)＝ae(t)yj(t)
4)对IMF分量中调频部分yj(t)进行Hilbert变换：
Figure GDA0002715474440000032
上式中，p为柯西主值。由此可得出IMF分量θ(t)的瞬时相位：
Figure GDA0002715474440000033
瞬时相位的导数即为IMF分量ci(t)的瞬时频率：
Figure GDA0002715474440000034
通过上述的方法，利用IMF分量的调频分量求取瞬时频率，利用其调幅分量求取瞬时振幅，弥补了常规Hilbert变换的不足，得到的瞬时属性能够更准确地反映信号的特征。但是在拟合经验包络信号时，由于IMF分量的绝对值|ci(t)|相比于|ci(t)|的局部极大值序列，具有较强的时变性，采用AR模型直接预测信号两端的数据会产生很大的误差。
在HHT的经验模式分解(EMD)过程中，存在模态混叠、端点效应、包络拟合误差以及虚假IMF分量问题，在现有技术中，解决模态混叠的方法存在计算复杂度高问题；在解决端点效应中信号两端存在较大的预测误差问题；计算均值曲线时，由于需要拟合上下两个包络，会产生较大的误差；在剔除虚假IMF分量问题采用单一指标筛选时，会产生筛选不准确的问题。而在第二步进行希尔伯特变换时，存在瞬时频率误差问题，根据上文阐述，归一化希尔伯特变换过程中做端点预测时存在较大的误差，所以需要做更精确的处理。
针对这些问题，本文提出解决这些问题的有效方法。
发明内容
本发明所要解决的技术问题是，提供一种能够克服现有HHT方法所存在的模态混叠、端点效应、包络拟合误差和瞬时频率误差问题，得到更准确的信号特征的改进的HHT方法。
本发明所采用的技术方案是：一种改进的HHT方法，包括如下步骤：
1)对原始信号挑选异常的极值点提取间歇信号，再运用自回归模型包络预测方法弥补异常的极值点，在波形上消除间歇信号抑制了模态混叠，从而得到修复后的信号，运用Sliding-fastBSpline-EMD方法分解信号，得到所有IMF分量；
2)根据得到的所有IMF分量与原始信号，计算每一个IMF分量的相关系数、互信息量和能量比值三组指标，并将每一个IMF分量的三组指标与三组指标相对应的设定阈值进行比较，根据比较结果对该IMF分量进行筛选，最后得到能够真实反映信号特征的一组IMF分量；
3)对所述一组IMF分量作基于自回归模型包络预测的归一化希尔伯特变换，根据希尔伯特变换得到IMF分量的瞬时相位，对瞬时相位求导，最终得到信号的瞬时频率。
步骤1)包括：
(1)挑选异常的极值点，根据非线性、非平稳信号中相邻极大值和极小值点符号是否相同，以及极值点是否为零来判断极值点是否异常，当相邻极大值和极小值符号同为正数时，极小值对应的点为异常极值点，当相邻极大值和极小值符号同为负数时，极大值对应的点为异常极值点；
(2)通过三次样条插值的方法，依据与异常极值点最近的3个以上连续的极值点，得到异常极值点的预测值，并用所述预测值替代原来的异常极值点，恢复一个相对正常的极值点序列，从而得到修复后的信号；
(3)运用Sliding-fastBSpline-EMD方法对修复后的信号进行EMD分解，得到所有IMF分量。
第(3)步所述的对修复后的信号进行EMD分解，包括：
(3.1)设修复后的信号的数据量长度为L，窗口长度为M，窗口每次滑动的长度为k，取M＝2d，d的取值范围为
Figure GDA0002715474440000041
以及滑动窗口数目为
Figure GDA0002715474440000042
但是由于在处理的信号中存在端点效应，取信号的前k(N-1)和后k(N-1)+mod(L,k)作为边界延拓数据；
其中，N为一个窗口长度中包含的滑动窗口个数。
(3.2)对于每一个滑动窗口，进行如下步骤：
(3.2a)设第i个滑动窗口在修复后的信号中的起点和终点分别是：start＝(i-1)·k+1和end＝start+M-1；
(3.2b)对滑动窗口内长度为M的数据采用fastBSpline-EMD方法分解，获得一组IMF分量：imfi,j，其中，j＝1,2,…,P，P为滑动窗口中信号分解得到的IMF分量的个数；
(3.2c)将每一个分量存入到一个P×kW的数组IMP中，公式如下：
IMPj,start～end＝IMPj,start～end+imfi,j；
式中，IMPj,start～end表示数组IMP的第j行，第start至end列的元素，imfi,j表示第j个IMF分量的值，数组IMP的每一行，即为1个IMF分量的值；
(3.2d)重复第(3.2a)～(3.2c)，直至最后一个滑动窗口；
(3.3)对存有所有滑动窗口信号的IMF分量累积数据的数组IMF求平均，得到所有IMF分量，公式如下：
Figure GDA0002715474440000043
式中，IMPj,k(N-1)～k(W-N+1)表示数组IMP的第j行，第k(N-1)至k(W-N+1)列的元素。
步骤2)所述的计算每一个IMF分量的互相关系数、互信息量和能量比值三组指标如下：
(1)互相关系数：
Figure GDA0002715474440000051
其中rj(t)是第j个IMF分量cj(t)与原始信号x(t)的互相关系数，P是信号经过EMD分解得到的所有IMF分量的个数。
(2)互信息量：
对于信号的第j个IMF分量cj(t)与原始信号x(t)
Figure GDA0002715474440000052
其中p(cj)是第j个IMF分量cj(t)的边缘概率分布，p(x)是原始信号x(t)的边缘概率分布，p(cj,x)是第j个IMF分量cj(t)与原始信号x(t)的联合概率分布。
(3)能量比值：
Figure GDA0002715474440000053
其中，N表示一个信号中采样点的个数，Enj是第j个IMF分量cj(t)的总能量，En原始信号x(t)的总能量，sj表示第j个IMF分量的总能量占原始信号总能量的比值，通过sj的值的大小来判断第j个IMF分量所占能量是否是虚假分量。
步骤2)所述的比较是：
设定互相关系数、互信息量和能量比值三组指标与对应的设定阈值比较结果分别是p1，p2和p3，当pz为1时，表示相应的指标大于该指标所对应的阈值，否则为0，其中z为1、2、3；
所述的筛选是：当p1+p2+p3≥2，则输出out＝1，否则out＝0，当输出为1，则表明所比较的IMF分量是真实分量，否则表明所比较的为虚假分量进行剔除。
步骤3)包括：
(1)对每一个待处理的IMF分量cj(t)求绝对值，记为|cj(t)|，找到|cj(t)|中所有的局部极大值；
(2)采用自回归模型对|cj(t)|进行包络预测，预测两端对应的极大值点，并将包含两端极大值点的|cj(t)|记为|cj′(t)|；
(3)使用三次样条插值方法对|cj′(t)|所有的局部极大值进行拟合，得到经验包络信号e1(t)；
(4)利用经验包络信号e1(t)对IMF分量cj(t)进行归一化处理，公式如下：
Figure GDA0002715474440000054
对y1(t)进行检验，若y1(t)≤1，则继续计算yh(t)＝yh-1(t)/ej(t)，其中h一次取值2,3,4，j＝h；直到计算结果满足yh(t)≤1的条件或者h＝4为止，此时，得到IMF分量cj(t)的调频部分为yh(t)；
(5)IMF分量cj(t)的调幅部分ae(t)用下式表示：
ae(t)＝e1(t)e2(t)…ej(t)
把IMF分量cj(t)用下式的形式表达：
cj(t)＝ae(t)yh(t)
(6)对IMF分量cj(t)中的调频部分yh(t)进行希尔伯特变换：
Figure GDA0002715474440000061
上式中，p为柯西主值，由此得出IMF分量cj(t)的瞬时相位θ(t)：
Figure GDA0002715474440000062
瞬时相位的导数即为IMF分量cj(t)的瞬时频率：
Figure GDA0002715474440000063
第(2)步包括：
若|cj(t)|的前q个数据{|cj(t-q)|,|cj(t-q+1)|,…,|cj(t-1)|}已知，根据所述的前q个数据对|cj(t)|的预测值记为
Figure GDA0002715474440000064
则有下列的自回归模型：
Figure GDA0002715474440000065
通过Levinson-Durbin算法求出AR模型中的系数{ak}，根据AR模型求出|cj(t)|的预测值，利用序列{)cj(t-q+1)|,…,)cj(t)|}对t+1时刻|cj(t+1)|的预测值
Figure GDA0002715474440000066
进行预测，依次类推求出以后任意时刻|cj(t)|的值；
利用自回归模型对|cj(t)|的两端极大值进行预测，构建新的IMF分量的绝对值|cj′(t)|。
本发明的一种改进的HHT方法，实现对非线性、非平稳的信号做更简单准确的时频分析，从而得到更准确的信号特征。本发明的方法能够有效地解决的抑制模态混叠、端点效应、包络拟合等问题带来的误差，同时减少Hilbert变换求取瞬时频率的误差问题。本发明的改进NHT相比于常规的NHT具有计算量小、准确的特点。
附图说明
图1是本发明一种改进的HHT方法的流程图；
图2是本发明中对信号进行Sliding-fastBSpline-EMD分解示意图。
具体实施方式
下面结合实施例和附图对本发明的一种改进的HHT方法做出详细说明。
如图1所示，本发明的一种改进的HHT方法，包括如下步骤：
1)对原始信号挑选异常的极值点提取间歇信号，再运用自回归模型包络预测方法弥补异常的极值点，在波形上消除间歇信号抑制了模态混叠，从而得到修复后的信号，运用Sliding-fastBSpline-EMD方法分解信号，得到所有IMF分量；如图2所示，包括：
(1)挑选异常的极值点，根据非线性、非平稳信号中相邻极大值和极小值点符号是否相同，以及极值点是否为零来判断极值点是否异常，当相邻极大值和极小值符号同为正数时，极小值对应的点为异常极值点，当相邻极大值和极小值符号同为负数时，极大值对应的点为异常极值点；
(2)通过三次样条插值的方法，依据与异常极值点最近的3个以上连续的极值点，得到异常极值点的预测值，并用所述预测值替代原来的异常极值点，恢复一个相对正常的极值点序列，从而得到修复后的信号；
(3)运用Sliding-fastBSpline-EMD方法对修复后的信号进行EMD分解，得到所有IMF分量。所述的对修复后的信号进行EMD分解，包括：
(3.1)设修复后的信号的数据量长度为L，窗口长度为M，窗口每次滑动的长度为k，取M＝2d，d的取值范围为
Figure GDA0002715474440000071
以及滑动窗口数目为
Figure GDA0002715474440000072
但是由于在处理的信号中存在端点效应，取信号的前k(N-1)和后k(N-1)+mod(L,k)作为边界延拓数据；
其中，N为一个窗口长度中包含的滑动窗口个数。
(3.2)对于每一个滑动窗口，进行如下步骤：
(3.2a)设第i个滑动窗口在修复后的信号中的起点和终点分别是：start＝(i-1)·k+1和end＝start+M-1；
(3.2b)对滑动窗口内长度为M的数据采用fastBSpline-EMD方法分解，获得一组IMF分量：imfi,j，其中，j＝1,2,…,P，P为滑动窗口中信号分解得到的IMF分量的个数；
(3.2c)将每一个分量存入到一个P×kW的数组IMP中，公式如下：
IMPj,start～end＝IMPj,start～end+imfi,j；
式中，IMPj,start～end表示数组IMP的第j行，第start至end列的元素，imfi,j表示第j个IMF分量的值，数组IMP的每一行，即为1个IMF分量的值；
(3.2d)重复第(3.2a)～(3.2c)，直至最后一个滑动窗口；
(3.3)对存有所有滑动窗口信号的IMF分量累积数据的数组IMF求平均，得到所有IMF分量，公式如下：
Figure GDA0002715474440000073
式中，IMPj,k(N-1)～k(W-N+1)表示数组IMP的第j行，第k(N-1)至k(W-N+1)列的元素。
2)根据得到的所有IMF分量与原始信号，计算每一个IMF分量的相关系数、互信息量和能量比值三组指标，并将每一个IMF分量的三组指标与三组指标相对应的设定阈值进行比较，根据比较结果对该IMF分量进行筛选，最后得到能够真实反映信号特征的一组IMF分量；
所述的计算每一个IMF分量的互相关系数、互信息量和能量比值三组指标如下：
(1)互相关系数：
Figure GDA0002715474440000074
其中rj(t)是第j个IMF分量cj(t)与原始信号x(t)的互相关系数，P是信号经过EMD分解得到的所有IMF分量的个数。
(2)互信息量：
对于信号的第j个IMF分量cj(t)与原始信号x(t)
Figure GDA0002715474440000081
其中p(cj)是第j个IMF分量cj(t)的边缘概率分布，p(x)是原始信号x(t)的边缘概率分布，p(cj,x)是第j个IMF分量cj(t)与原始信号x(t)的联合概率分布。
(3)能量比值：
Figure GDA0002715474440000082
其中，N表示一个信号中采样点的个数，Enj是第j个IMF分量cj(t)的总能量，En原始信号x(t)的总能量，sj表示第j个IMF分量的总能量占原始信号总能量的比值，通过sj的值的大小来判断第j个IMF分量所占能量是否是虚假分量。
所述的比较是：
设定互相关系数、互信息量和能量比值三组指标与对应的设定阈值比较结果分别是p1，p2和p3，当pz为1时，表示相应的指标大于该指标所对应的阈值，否则为0，其中z为1、2、3；
所述的筛选：当p1+p2+p3≥2，则输出out＝1，否则out＝0，当输出为1，则表明所比较的IMF分量是真实分量，否则表明所比较的为虚假分量进行剔除。
3)对所述一组IMF分量作基于AR模型包络预测的归一化希尔伯特变换(NHT)，根据希尔伯特变换得到IMF分量的瞬时相位，对瞬时相位求导，最终得到信号的瞬时频率。
改进的NHT变换，使用自回归模型包络预测方法代替原来提到的自回归模型的信号预测延拓方法。利用经验包络信号根据归一化检验规则对IMF分量进行归一化处理，然后对IMF分量中的调频部分进行Hilbert变换，最终经过求导得到各个分量对应的瞬时频率。
该步骤包括：
(1)对每一个待处理的IMF分量cj(t)求绝对值，记为|cj(t)|，找到|cj(t)|中所有的局部极大值；
(2)采用自回归模型对|cj(t)|进行包络预测，预测两端对应的极大值点，并将包含两端极大值点的|cj(t)|记为|cj′(t)|；包括：
若|cj(t)|的前q个数据{|cj(t-q)|,|cj(t-q+1)|,…,|cj(t-1)|}已知，根据所述的前q个数据对|cj(t)|的预测值记为
Figure GDA0002715474440000083
则有下列的自回归模型：
Figure GDA0002715474440000084
通过Levinson-Durbin算法求出自回归模型中的系数{ak}，根据自回归模型求出|cj(t)|的预测值，利用序列{|cj(t-q+1)|,…,|cj(t)|}对t+1时刻|cj(t+1)|的预测值
Figure GDA0002715474440000085
进行预测，依次类推求出以后任意时刻)cj(t)|的值；
利用上述自回归模型对|cj(t)|的两端极大值进行预测，构建新的IMF分量的绝对值|cj′(t)|。
(3)使用三次样条插值方法对|cj′(t)|所有的局部极大值进行拟合，得到经验包络信号e1(t)；
(4)利用经验包络信号e1(t)对IMF分量cj(t)进行归一化处理，公式如下：
Figure GDA0002715474440000091
对y1(t)进行检验，若y1(t)≤1，则继续计算yh(t)＝yh-1(t)/ej(t)，其中h一次取值2,3,4，j＝h；直到计算结果满足yh(t)≤1的条件或者h＝4为止，此时，得到IMF分量cj(t)的调频部分为yh(t)；
(5)根据第(4)步，cj(t)的调幅部分ae(t)用下式表示：
ae(t)＝e1(t)e2(t)…ej(t)
把IMF分量cj(t)用下式的形式表达：
cj(t)＝ae(t)yh(t)
(6)对IMF分量中调频部分yh(t)进行希尔伯特变换：
Figure GDA0002715474440000092
上式中，p为柯西主值，由此得出IMF分量的瞬时相位θ(t)：
Figure GDA0002715474440000093
瞬时相位的导数即为IMF分量cj(t)的瞬时频率：
Figure GDA0002715474440000094
。

Claims (7)

1.一种改进的HHT方法，应用于信号处理，其特征在于，包括如下步骤：
1)对原始信号挑选异常的极值点提取间歇信号，再运用自回归模型包络预测方法弥补异常的极值点，在波形上消除间歇信号抑制了模态混叠，从而得到修复后的信号，运用Sliding-fastBSpline-EMD方法分解信号，得到所有IMF分量；
2)根据得到的所有IMF分量与原始信号，计算每一个IMF分量的相关系数、互信息量和能量比值三组指标，并将每一个IMF分量的三组指标与三组指标相对应的设定阈值进行比较，根据比较结果对该IMF分量进行筛选，最后得到能够真实反映信号特征的一组IMF分量；
3)对所述一组IMF分量作基于自回归模型包络预测的归一化希尔伯特变换，根据希尔伯特变换得到IMF分量的瞬时相位，对瞬时相位求导，最终得到信号的瞬时频率。
2.根据权利要求1所述的一种改进的HHT方法，其特征在于，步骤1)包括：
(1)挑选异常的极值点，根据非线性、非平稳信号中相邻极大值和极小值点符号是否相同，以及极值点是否为零来判断极值点是否异常，当相邻极大值和极小值符号同为正数时，极小值对应的点为异常极值点，当相邻极大值和极小值符号同为负数时，极大值对应的点为异常极值点；
(2)通过三次样条插值的方法，依据与异常极值点最近的3个以上连续的极值点，得到异常极值点的预测值，并用所述预测值替代原来的异常极值点，恢复一个相对正常的极值点序列，从而得到修复后的信号；
(3)运用Sliding-fastBSpline-EMD方法对修复后的信号进行EMD分解，得到所有IMF分量。
3.根据权利要求2所述的一种改进的HHT方法，其特征在于，第(3)步所述的对修复后的信号进行EMD分解，包括：
(3.1)设修复后的信号的数据量长度为L，窗口长度为M，窗口每次滑动的长度为k，取M＝2d，d的取值范围为
Figure FDA0002743582210000011
以及滑动窗口数目为
Figure FDA0002743582210000012
但是由于在处理的信号中存在端点效应，取信号的前k×(N-1)和后k×(N-1)+mod(L,k)作为边界延拓数据；
其中，N为一个窗口长度中包含的滑动窗口个数；
(3.2)对于每一个滑动窗口，进行如下步骤：
(3.2a)设第i个滑动窗口在修复后的信号中的起点和终点分别是：start＝(i-1)×k+1和end＝start+M-1；
(3.2b)对滑动窗口内长度为M的数据采用fastBSpline-EMD方法分解，获得一组IMF分量：imfi,j，其中，j＝1,2,…，P，P为滑动窗口中信号分解得到的IMF分量的个数；
(3.2c)将每一个分量存入到一个P×kW的数组IMP中，公式如下：
IMPj,start～end＝IMPj,start～end+imfi,j；
式中，IMPj,start～end表示数组IMP的第j行，第start至end列的元素，imfi,j表示第j个IMF分量的值，数组IMP的每一行，即为1个IMF分量的值；
(3.2d)重复第(3.2a)～(3.2c)，直至最后一个滑动窗口；
(3.3)对存有所有滑动窗口信号的IMF分量累积数据的数组IMF求平均，得到所有IMF分量，公式如下：
Figure FDA0002743582210000021
式中，IMPj,k×(N-1)～k×(W-N+1)表示数组IMP的第j行，第k×(N-1)至k×(W-N+1)列的元素。
4.根据权利要求1所述的一种改进的HHT方法，其特征在于，步骤2)所述的计算每一个IMF分量的互相关系数、互信息量和能量比值三组指标如下：
(1)互相关系数：
Figure FDA0002743582210000022
其中rj(t)是第j个IMF分量cj(t)与原始信号x(t)的互相关系数，P是信号经过EMD分解得到的所有IMF分量的个数，
Figure FDA0002743582210000023
是第j个IMF分量cj(t)的平均值，
Figure FDA0002743582210000024
是原始信号x(t)的平均值；
(2)互信息量：
对于信号的第j个IMF分量cj(t)与原始信号x(t)
Figure FDA0002743582210000025
其中p(cj)是第j个IMF分量cj(t)的边缘概率分布，p(x)是原始信号x(t)的边缘概率分布，p(cj,x)是第j个IMF分量cj(t)与原始信号x(t)的联合概率分布；
(3)能量比值：
Figure FDA0002743582210000026
其中，N表示一个信号中采样点的个数，Enj是第j个IMF分量cj(t)的总能量，En原始信号x(t)的总能量，sj表示第j个IMF分量的总能量占原始信号总能量的比值，通过sj的值的大小来判断第j个IMF分量所占能量是否是虚假分量。
5.根据权利要求1所述的一种改进的HHT方法，其特征在于，步骤2)所述的比较是：
设定互相关系数、互信息量和能量比值三组指标与对应的设定阈值比较结果分别是p1，p2和p3，当pz为1时，表示相应的指标大于该指标所对应的阈值，否则为0，其中z为1、2、3；
所述的筛选是：当p1+p2+p3≥2，则输出out＝1，否则out＝0，当输出为1，则表明所比较的IMF分量是真实分量，否则表明所比较的为虚假分量进行剔除。
6.根据权利要求1所述的一种改进的HHT方法，其特征在于，步骤3)包括：
(1)对每一个待处理的IMF分量cj(t)求绝对值，记为|cj(t)|，找到|cj(t)|中所有的局部极大值；
(2)采用自回归模型对|cj(t)|进行包络预测，预测两端对应的极大值点，并将包含两端极大值点的|cj(t)|记为|cj′(t)|；
(3)使用三次样条插值方法对|cj′(t)|所有的局部极大值进行拟合，得到经验包络信号e1(t)；
(4)利用经验包络信号e1(t)对IMF分量cj(t)进行归一化处理，公式如下：
Figure FDA0002743582210000031
对y1(t)进行检验，若y1(t)≤1，则结束归一化过程，直接得到IMF分量cj(t)的调频部分为y1(t)；若y1(t)>1，则继续计算yh(t)＝yh-1(t)/ej(t)，其中h依次取值2,3,4，j＝h；直到计算结果满足yh(t)≤1的条件或者h＝4为止，此时，得到IMF分量cj(t)的调频部分为yh(t)；
(5)IMF分量cj(t)的调幅部分ae(t)用下式表示：
ae(t)＝e1(t)e2(t)…ej(t)
把IMF分量cj(t)用下式的形式表达：
cj(t)＝ae(t)yh(t)
(6)对IMF分量cj(t)中的调频部分yh(t)进行希尔伯特变换：
Figure FDA0002743582210000032
上式中，p为柯西主值，由此得出IMF分量cj(t)的瞬时相位θ(t)：
Figure FDA0002743582210000033
瞬时相位的导数即为IMF分量cj(t)的瞬时频率：
Figure FDA0002743582210000034
7.根据权利要求6所述的一种改进的HHT方法，其特征在于，第(2)步包括：
若|cj(t)|的前q个数据{|cj(t-q)|,|cj(t-q+1)|,…,|cj(t-1)|}已知，根据所述的前q个数据对|cj(t)|的预测值记为
Figure FDA0002743582210000035
则有下列的自回归模型：
Figure FDA0002743582210000036
通过Levinson-Durbin算法求出AR模型中的系数ak，根据自回归模型求出|cj(t)|的预测值，利用序列{|cj(t-q+1)|,…,|cj(t)|}对t+1时刻|cj(t+1)|的预测值
Figure FDA0002743582210000037
进行预测，依次类推求出以后任意时刻|cj(t)|的值；
利用自回归模型对|cj(t)|的两端极大值进行预测，构建新的IMF分量的绝对值|cj′(t)|。
CN201711081152.9A 2017-11-06 2017-11-06 一种改进的hht方法 Active CN108009122B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711081152.9A CN108009122B (zh) 2017-11-06 2017-11-06 一种改进的hht方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711081152.9A CN108009122B (zh) 2017-11-06 2017-11-06 一种改进的hht方法

Publications (2)

Publication Number Publication Date
CN108009122A CN108009122A (zh) 2018-05-08
CN108009122B true CN108009122B (zh) 2021-01-05

Family

ID=62051319

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711081152.9A Active CN108009122B (zh) 2017-11-06 2017-11-06 一种改进的hht方法

Country Status (1)

Country Link
CN (1) CN108009122B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108957175B (zh) * 2018-06-15 2020-09-25 西安理工大学 基于改进的hht算法的电能质量扰动识别方法
CN109460614B (zh) * 2018-11-12 2023-06-30 广西交通科学研究院有限公司 基于瞬时带宽的信号时间-频率分解方法
CN109558863B (zh) * 2019-01-07 2020-11-06 哈尔滨工业大学（深圳） 脑电信号特征表示与提取方法、装置及存储介质
CN110514921B (zh) * 2019-07-22 2022-07-26 华南理工大学 一种电力电子变换器非平稳信号中非线性现象的识别方法
CN112270209A (zh) * 2020-09-29 2021-01-26 河南科技大学 一种联合收割机脱粒滚筒振动信号特征提取方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5983162A (en) * 1996-08-12 1999-11-09 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Computer implemented empirical mode decomposition method, apparatus and article of manufacture
US7747401B2 (en) * 2007-08-03 2010-06-29 Oracle International Corporation Fast intrinsic mode decomposition of time series data with sawtooth transform
US7941298B2 (en) * 2006-09-07 2011-05-10 DynaDx Corporation Noise-assisted data analysis method, system and program product therefor
CN104680011A (zh) * 2015-02-16 2015-06-03 燕山大学 基于amd的消除经验模态分解中模态混叠方法
CN105760347A (zh) * 2016-02-04 2016-07-13 福建工程学院 一种基于数据/极值联合对称延拓的hht端点效应抑制方法
CN106844293A (zh) * 2017-02-03 2017-06-13 中国铁道科学研究院 一种经验模态分解中模态混叠问题的自适应解耦方法
CN106940688A (zh) * 2017-01-08 2017-07-11 广东工业大学 一种emd虚假分量识别的数学模型建模方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TWI460416B (zh) * 2011-03-28 2014-11-11 Univ Nat Taiwan 機械系統狀態之判斷方法及判斷裝置

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5983162A (en) * 1996-08-12 1999-11-09 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Computer implemented empirical mode decomposition method, apparatus and article of manufacture
US7941298B2 (en) * 2006-09-07 2011-05-10 DynaDx Corporation Noise-assisted data analysis method, system and program product therefor
US7747401B2 (en) * 2007-08-03 2010-06-29 Oracle International Corporation Fast intrinsic mode decomposition of time series data with sawtooth transform
CN104680011A (zh) * 2015-02-16 2015-06-03 燕山大学 基于amd的消除经验模态分解中模态混叠方法
CN105760347A (zh) * 2016-02-04 2016-07-13 福建工程学院 一种基于数据/极值联合对称延拓的hht端点效应抑制方法
CN106940688A (zh) * 2017-01-08 2017-07-11 广东工业大学 一种emd虚假分量识别的数学模型建模方法
CN106844293A (zh) * 2017-02-03 2017-06-13 中国铁道科学研究院 一种经验模态分解中模态混叠问题的自适应解耦方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A NEW VIEW OF NONLINEAR;Norden E. Huang等;《arjournals.annualreviews.org》;19990331;第417-457页 *
DETECTION OF THE VORTEX SIGNAL IN VIBRATION ENVIRONMENT;HONG-JUN SUN等;《Proceedings of the 2010 International Conference on Wavelet Analysis and Pattern Recognition》;20100711;第216-221页 *
EMD间歇信号的检测和提取方法;胡重庆等;《数据采集与处理》;20080131;第23卷(第1期);第108-111页 *
Hilbert‐Huang变换中的模态混叠问题;曹 莹等;《振动、测试与诊断》;20160630;第36卷(第3期);第518-523页 *

Also Published As

Publication number Publication date
CN108009122A (zh) 2018-05-08

Similar Documents

Publication Publication Date Title
CN108009122B (zh) 一种改进的hht方法
US7401008B2 (en) Method, computer program, and system for intrinsic timescale decomposition, filtering, and automated analysis of signals of arbitrary origin or timescale
CN109034127B (zh) 一种频谱异常检测方法、装置和电子设备
CN109727446A (zh) 一种用电数据异常值的识别与处理方法
CN109239763B (zh) 模拟核衰变过程的模拟能谱曲线仿真方法
CN117434153B (zh) 基于超声波技术的道路无损检测方法及系统
CN113392732A (zh) 一种局部放电超声信号抗干扰方法及系统
CN107561420A (zh) 一种基于经验模态分解的电缆局部放电信号特征向量提取方法
CN108415880B (zh) 一种基于样本熵与小波变换的线损特性分析方法
CN109065176B (zh) 一种血糖预测方法、装置、终端和存储介质
WO2015074884A1 (en) Methods and systems for wavelet based representation
CN110751400B (zh) 一种风险评估方法及装置
CN112151067A (zh) 一种基于卷积神经网络的数字音频篡改被动检测方法
CN116979970A (zh) 一种漏磁数据的压缩和重构方法、系统、电子设备及介质
CN116758922A (zh) 一种用于变压器的声纹监测与诊断方法
CN116223627A (zh) 一种钢轨探伤数据降噪装置
CN113598784B (zh) 心律失常检测方法及系统
CN115859048A (zh) 一种局放信号的噪声处理方法及装置
JPH0573093A (ja) 信号特徴点の抽出方法
CN112699614A (zh) 基于XGBoost的序列预测模型构建、降水趋势预测方法及装置
Maulana The utilization of Gaussian filter method on voice record frequency noise
CN112149833A (zh) 基于机器学习的预测方法、装置、设备和存储介质
CN113598810B (zh) 一种基于分割网络的胎心率基线自动计算方法
CN115270906A (zh) 基于电网频率深浅层特征融合的数字音频篡改被动检测方法及装置
Hanke An ϵ-free a posteriori stopping rule for certain iterative regularization methods

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant