CN103675444A - 一种高精度的时频分析方法 - Google Patents

Abstract

本发明提供了一种高精度的时频分析方法，属于数字信号处理领域。所述方法包括以下步骤：(1)输入原始数字信号x(t)，然后对所述原始数字信号x(t)进行EMD分解，将原始数字信号x(t)分解成IMF分量的集合，即

其中c_i(t)为IMF分量；(2)分别对每个IMF分量c_i(t)进行MP分解，得到每个IMF分量的时频谱；(3)将步骤(2)得到的各个IMF分量的时频谱相互叠加，得到原始数据信号x(t)的高时频分辨率的时频谱。本发明方法避免了直接对复合信号进行MP分解匹配时频原子的不准确性，通过对各IMF分量进行MP分解能够得到准确性更高、时频分辨率更高、信息更丰富的时频谱。

Description

一种高精度的时频分析方法

技术领域

本发明属于数字信号处理领域，具体涉及一种针对数字信号的高精度的时频分析方法。

背景技术

时频分析是分析非平稳信号的传统方法，它能将一维的时域或频域信号变换到二维时频平面上，得到信号的时频谱，可以在时频域中提取并区分信号分量，能够更清楚的了解信号的特征。时频分析方法一般分为线性时频分析方法、双线性时频分析方法和参数化时频分析方法等。典型的线性时频分析方法包含短时傅里叶变换、小波变换、S变换等；由于线性时频分析方法受Heisenberg测不准原理的约束，信号的时间分辨率和频率分辨率为一对矛盾体，不能同时达到较高的时间分辨率和较高的频率分辨率。典型的双线性时频分析方法为Wigner-Ville分布(即WVD)；由于WVD中不含任何窗函数，避免了线性时频分析中时间分辨率和频率分辨率矛盾的问题；在分解单分量平稳信号时具有很高的时频分辨率；但是WVD为双线性的，在分解多分量非平稳信号时会出现严重的交叉项。

最典型的参数化时频分析方法为匹配追踪法。匹配追踪法为S.Mallat与Z.Zhang于1993年提出的，它是通过对时频原子进行伸缩、时移、频移和相移得到一个原子库，同时根据最大匹配投影原理寻找最佳原子的线性组合，以达到对原始信号进行分解，具有较高的时频分辨率。但是MP算法是一种贪婪算法，它是通过不断地迭代寻找局部最佳原子的方法来实现对原信号的不断逼近过程。为了提高计算效率，YangHua Wang提出了在迭代过程中全局粗粒度预测和局部最优相结合的方法，充分考虑了尺度的重要性，提出了尺度滤波的思想，提高了算法的计算效率。

但是对于复合信号，由于波与波相互叠加形成了复合波，直接对复合波进行匹配搜索往往得不到单一的时频原子。如在某一个时间段，彼此间时移较小的几个单波相互叠加。若直接利用MP算法进行匹配，匹配的时频原子便为一个复合的时频原子，为几个单波相叠加共同作用的波；按照匹配的思想需要继续迭代，需要从原始信号中减去搜索的时频原子，那么剩下的波形便不是一个完整的波形，而是复合波的剩余波，它包含每一个波的部分信息。那么随着不断的搜索，匹配的时频原子越来越失真、不具任何意义。同时，直接利用对非线性的复合信号经希尔伯特变换计算的属性来预测时频原子的时间、频率和相位是没有意义的，不能准确的定位单组份信号的时间和频率等相关信息，导致最终计算的时频谱不准确。

瞬时频率为解析信号瞬时相位的导数，为了得到有意义的瞬时频率，E.N.Huang等人对瞬时频率进行了充分的分析，认为只有本征模函数经希尔伯特变换计算的瞬时频率才有意义。在此基础上，提出了一种全新的信号时频分析方法——希尔伯特黄变换(Hilbert-Huang Transform，即HHT)。HHT最关键的技术为经验模态分解(Empirical Mode Decomposition，即EMD)，即对时间序列进行经验模态分解，分解成本征模函数组。经验模态分解从本质上讲是对信号进行平稳化的过程，将信号中客观存在的、不同尺度的波动分解出来，得到一系列不同尺度数据序列，即本征模函数。

发明内容

本发明的目的在于解决目前匹配追踪时频分解中时频原子搜索的准确性问题，提供一种高精度的时频分析方法，充分考虑了IMF(Intrinsic ModeFunction)分量所特有的性质，以及希尔伯特变换对实信号的要求；采取的是对IMF分量进行MP分解(Matching Pursuit Decomposition，即MP分解)，能够更准确地定位时频原子的位置、频率、相位及振幅等信息，能够检测出更丰富的频率成分。采用时频分辨率较高的WVD计算时频原子的时频谱，能够得到高时频分辨率的时频谱。

本发明是通过以下技术方案实现的：

所述方法包括以下步骤：

(1)输入原始数字信号x(t)，然后对所述原始数字信号x(t)进行EMD分解，将原始数字信号x(t)分解成IMF分量的集合，即

其中c_i(t)为IMF分量；

(2)分别对每个IMF分量c_i(t)进行MP分解，得到每个IMF分量的时频谱；

(3)将步骤(2)得到的各个IMF分量的时频谱相互叠加(即同一时刻、同一频率上的值相加)，得到原始数据信号x(t)的高时频分辨率的时频谱。

所述步骤(1)具体包括以下步骤：

(11)初始化，r₀(t)＝x(t)，令i＝1；

(12)筛选第i个IMF分量，具体包括以下步骤：

(12a)初始化，h₀(t)＝r_i-1(t)，j＝1；

(12b)根据信号h_j-1(t)的波形，找出h_j-1(t)的局部极大值点和局部极小值点的位置和振幅值，并采用边界处理方法分别将所述局部极大值点和局部极小值点向两边延拓1-2个极值点；

(12c)利用差值方法对h_j-1(t)的局部极小值点进行差值，得到h_j-1(t)的下包络线；对h_j-1(t)的局部极大值点进行差值，得到h_j-1(t)的上包络线；

(12d)计算步骤(12c)得到的上包络线和下包络线的平均值m_j-1(t)；

(12e)做减法：h_j(t)＝h_j-1(t)-m_j-1(t)；

(12f)如果满足终止条件，则令c_i(t)＝h_j(t)，然后转入步骤(13)；否则，j＝j+1，转到(12b)；

所述边界处理方法可以采用镜像延拓法、神经网络预测法或波形匹配法等。针对不同类型的信号采取不同的处理方法能达到不同的效果。

对局部极小值点和局部极大值点进行差值可以采用三次样条差值、B样条函数的线性组合法、局部滤波法、分段幂函数插值法或高阶样条拟合法等方法。

所述终止条件可以采用满足仿柯两准则或满足极值点数目与零点数目相等或最多相差一个等。

(13)求剩余信号：r_i(t)＝r_i-1(t)-c_i(t)；

(14)判断r_i(t)的极值点的个数是否仍然超过两个，如果是，则i＝i+1，转到步骤(12)，如果否，则分解结束，r_i(t)即为剩余分量。

所述步骤(2)是这样实现的：首先进行MP分解，即搜索组成一个IMF分量的时频原子，然后计算每个时频原子的时频谱，最后将所有时频原子的时频谱相叠加便得到该IMF分量的时频谱，具体包括以下步骤：

(21)初始化，y_j(t)＝c_i(t)，令j＝1；

(22)对y_j(t)进行希尔伯特变换，生成复信号，然后分别求解复信号的瞬时振幅、瞬时频率和瞬时相位；其中u_n为瞬时振幅最大值所对应的时间；f_n为该时间对应的瞬时频率；φ_n为该时间对应的瞬时相位；

对数字信号进行希尔伯特变换、计算瞬时振幅、瞬时相位和瞬时频率都是公知的常用方法，在三瞬分析中经常使用，其变换和计算方法在公知教科书中都有详细说明。

在一组固定的u_n、f_n和φ_n情况下，通过计算最优化公式得到尺度σ_n；其中D＝{g_r(t)}_r∈Γ为时频原子字典，

是函数R⁽ⁿ⁾s和

的内积，

通过这一步便得到第j个时频原子四个参数u_n、f_n、φ_n和σ_n的初始值；

(23)利用所述最优化公式在局域内寻找r_n＝{u_n，σ_n，f_n，φ_n}的最优值，即公式

达到最大时，对应的一组参数便是最优的r_n；搜索范围为[r_n-Δr，r_n+Δr]，其中Δr＝(Δu，Δσ，Δf，Δφ)，Δu为时间偏移量、Δσ为尺度偏移量、Δf为频率偏移量、Δφ为相位偏移量；

(24)在上面两步的基础上，便可以预测到最优的时频原子

的振幅为a_n，即

然后计算剩余信号

令j＝j+1，重复迭代步骤(22)到(24)N次，直到剩余量y_n(t)的最大能量小于原始信号最大能量的10％或更低，得到该IMF分量的N个时频原子，y_n(t)可视为IMF分量c_i(t)的噪音部分；

(25)利用WVD分别计算所述N个时频原子的时频谱，然后将所述N个时频原子的时频谱叠加起来，便得到该IMF分量的时频谱；

(26)对每个IMF分量进行步骤(21)至(25)的处理得到每个IMF分量的时频谱。

与现有技术相比，本发明的有益效果是：

(1)本发明方法是一种高分辨率的时频分析方法；

(2)本发明方法利用EMD分解将复合信号分解成单分量的窄带信号，而希尔伯特变换对于这种单分量信号才有意义；

(3)本发明方法对单分量信号进行MP分解能够得到更精确的时频原子和更丰富的频率信息，避免了相邻波叠加造成时频原子匹配的不准确性问题，能够得到时频分辨率更高的时频谱。

附图说明

图1是本发明高精度的时频分析方法的步骤框图。

图2-1是利用一个频率跃变信号测试瞬时频率准确性的实施例中的各分量及合成信号。

图2-2是利用一个频率跃变信号测试瞬时频率准确性的实施例中的合成信号的瞬时频率。

图2-3是利用一个频率跃变信号测试瞬时频率准确性的实施例中的各分量的瞬时频率。

图3是一道合成地震信号及其几个分量示意图。

图4-1是利用一道合成地震信号测试瞬时频率准确性的实施例中的各分量及合成信号。

图4-2是利用一道合成地震信号测试瞬时频率准确性的实施例中的合成信号的瞬时频率。

图4-3是利用一道合成地震信号测试瞬时频率准确性的实施例中的各分量的瞬时频率。

图5是利用MP算法对合成地震信号进行匹配搜索的时频原子示意图。

图6-1是直接对组成信号的各时频原子计算WVD叠加起来的时频谱。

图6-2是直接对合成信号进行MP分解计算的时频谱。

图6-3是利用本发明方法计算的时频谱。

图7是一道实际的地震资料及其分量示意图。

图8-1是实际地震资料及其分量的瞬时频率对比实施例中的原始信号的瞬时频率。

图8-2是实际地震资料及其分量的瞬时频率对比实施例中的各分量的瞬时频率。

图9-1是直接对地震资料进行MP分解得到的时频谱图。

图9-2是利用本发明方法计算得到的时频谱图。

具体实施方式

下面结合附图对本发明作进一步详细描述。本发明的范围不受这些实施的限制，本发明的范围在权利要求书中已提出。

一种高精度时频分析方法，所述方法首先对一道时间序列数字信号(记为x(t))利用EMD方法分解成本征模函数(IMF)的集合；然后对每一个IMF分量进行匹配追踪(MP)分解，对匹配的时频原子进行WVD分析，得到时频原子的时频谱，将所有时频原子的时频谱叠加起来得到每个IMF分量的时频谱；最后将所有IMF分量的时频谱叠加起来便得到原始数字信号的高精度时频谱。本发明方法避免了直接对复合地震信号进行MP分解匹配时频原子的不准确性，通过对各IMF分量进行MP分解能够得到准确性更高、时频分辨率更高、信息更丰富的时频谱。

如图1所示，本发明高精度的时频分析方法包括以下步骤：

第一步：输入一道时间序列的原始数字信号x(t)，对原始数字信号x(t)进行EMD分解，分解成若干个IMF的集合，即

c_i(t)即为IMF分量。

第二步：对每一个IMF分量c_i(t)进行MP分解，分解成一组时频原子的集合，利用WVD计算各时频原子的时频谱，将所有时频原子的时频谱叠加起来，便得到每一个IMF分量的时频谱。

第三步：将各个IMF分量的时频谱叠加起来，得到最终的高时频分辨率的时频谱。

下面通过几个实施例来说明本发明方法的效果：

1，利用频率跃变信号测试瞬时频率准确性的实施例如下：

图2-1的三条曲线从上到下分别为信号1、信号2和由信号1+信号2组成的频率跃变信号(下称合成信号)；图2-2为直接对合成信号经希尔伯特变换计算的瞬时频率曲线；图2-3中的虚线为信号1的瞬时频率曲线，实线为信号2的瞬时频率曲线。从叠加性原理可知，合成信号由两个信号组成，两个分量在同一时刻具有不同的频率，那么合成信号在每一点也应该有两个不同频率的分量。但是从图2-2可以看出，直接对合成信号经希尔伯特变换计算的瞬时频率不能全部展现出信号中所包含的频率信息，表现的是一种混合信息。导致这种结果的原因在于合成信号不是单分量信号，不满足希尔伯特变换的要求，计算的瞬时频率不具有实际的意义，也是不准确的。

如果直接对合成信号进行MP分解，搜索的时频原子也是一种复合的时频原子，时频原子的时移、频率、相位、能量大小等参数都是不准确的，导致时频谱也不准确。而通过EMD分解能够将合成信号分解成独立的两个单分量信号，分解的两个单分量信号与合成的信号1和信号2一致，如果分别对两个单分量信号进行MP分解，能够匹配更准确的时频原子，得到精度更高的时频谱。

2，利用一道合成地震信号测试瞬时频率准确性的实施例如下：

图3显示了五个地震信号，从左到右分别为信号1、信号2、信号3、信号4和信号5，其中信号5为信号1、信号2、信号3和信号4叠加的结果，简称合成信号。信号1、信号2、信号3和信号4都为单分量信号，合成信号为复合信号。

图4-1与图3一致，图4-2显示的是直接对合成信号计算的瞬时频率。图4-3表示的是分别对四个单分量信号计算的瞬时频率，其中虚线表示的是信号1的瞬时频率，点线表示的是信号2的瞬时频率，实线表示的是信号3的瞬时频率，加线表示的是信号4的瞬时频率。对比图4-2和图4-3可以看出，在没有复合波的地方直接对合成信号计算瞬时频率与对应的分量计算的瞬时频率一致；而在时频原子相互叠加的部分，直接对合成信号计算瞬时频率不能完全展现出信号中包含的频率成分，因为合成信号在该时间段为多分量信号。

图5所示的是直接对合成信号做MP分解匹配的时频原子，与图3所示的合成信号的各分量相比，0-100ms之间只包含一个时频原子，匹配较准确；200-450ms之间包含两个频率相差较大、时移一致的时频原子，匹配也较准确。而在500-600ms之间合成信号应该包含两个频率、相位相同、时移相差较小的分量；但是从图5中可以看出，在这个时间段搜索的时频原子包含第7、第8和第13个时频原子，第7和第8个时频原子相位、能量不准确，而且还搜出信号中不存在的第13个时频原子。而在时移和频率都比较接近的时间范围内，如在650-850ms与850-1000ms之间只能搜索出复合的时频原子(如第2、3、4个时频原子)和一些原始信号中不存在的时频原子(如第9、10、11、12、14、15、16个时频原子)。

对合成地震信号进行MP分解得到的时频谱对比图如图6-1、图6-2和图6-3所示。其中图6-1表示的是直接对组成信号的各时频原子计算WVD叠加起来的时频谱，理论上讲这样的时频分辨率是最高的，不存在复合波；图6-2表示的是直接对合成信号进行MP分解计算的时频谱；图6-3为利用本发明计算的时频谱。对比图6-2和图6-3可知，利用本发明方法匹配的时频原子更准确，计算的时频谱分辨率更高，几平接近理论时频谱的分辨率。而直接对合成信号计算的时频谱时频分辨率不高，而且出现了很多原始信号中不存在的频率成分。

图7为一道实际的地震信号及其分量，分量为EMD对原始信号进行分解的结果。从下到上分别为地震信号、信号的IMF1、信号的IMF2、信号的IMF3以及信号的IMF4。从图中可以看出，原始地震信号不是一个单分量的信号，它包含了几个单分量信号。

实际地震信号的瞬时频率对比图如图8-1和图8-2所示。图8-1为直接对地震信号经希尔伯特变换计算的瞬时频率；图8-2为各分量的瞬时频率，其中虚线为IMF1的瞬时频率、点线为IMF2的瞬时频率、实线为IMF3的瞬时频率、点划线为IMF4的瞬时频率。对比图8-1和图8-2可知，EMD能够将复合波分解成单分量信号，有利于检测出更丰富的频率信息。

图9-1显示了直接对地震信号进行MP分解得到的时频谱，图9-2显示了利用本发明方法计算的时频谱。对比两图可知，利用本发明计算的时频谱能够识别出同一时刻或相邻时刻包含的多个不同频率的信号，分解出更丰富的分量信息。存在复合波的地方，利用EMD将其分解，得到单分量信号，消除了复合波之间的干扰，有利于匹配追踪匹配更准确地时频原子。如在1000ms-1500ms之间，直接对原始地震信号进行MP分解只能检测出约6个分量的信号，而利用本发明首先将复合信号进行分离，然后进行匹配，大约能检测出11个分量的信号，而且分量彼此间相邻较近。而利用传统的短时傅里叶变换、小波变换、S变换以及WVD等都不能得到这样的时频分辨率。

总之，本发明方法通过将复合信号通过EMD分解成单分量信号的集合，分别对各单分量信号进行MP分解，得到单分量信号的时频谱，最终将所有单分量信号的时频谱叠加起来得到高精度的时频谱。避免了直接对原始复合信号进行MP分解不能准确检测时频原子、时频谱计算不准确以及时频分辨率不高的问题。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (3)

1.一种高精度的时频分析方法，其特征在于：所述方法包括以下步骤：

(1)输入原始数字信号x(t)，然后对所述原始数字信号x(t)进行EMD分解，

将原始数字信号x(t)分解成IMF分量的集合，即

其中c_i(t)为IMF分量；

(2)分别对每个IMF分量c_i(t)进行MP分解，得到每个IMF分量的时频谱；

(3)将步骤(2)得到的各个IMF分量的时频谱相互叠加，得到原始数据信号x(t)的高时频分辨率的时频谱。

2.根据权利要求1所述的高精度的时频分析方法，其特征在于：所述步骤(1)具体包括以下步骤：

(11)初始化，r₀(t)＝x(t)，令i＝1；

(12)筛选第i个IMF分量，具体包括以下步骤：

(12a)初始化，h₀(t)＝r_i-1(t)，j＝1；

(12b)根据信号h_j-1(t)的波形，找出h_j-1(t)的局部极大值点和局部极小值点的位置和振幅值，并采用边界处理方法分别将所述局部极大值点和局部极小值点向两边延拓1-2个极值点；

(12c)利用差值方法对h_j-1(t)的局部极小值点进行差值，得到h_j-1(t)的下包络线；对h_j-1(t)的局部极大值点进行差值，得到h_j-1(t)的上包络线；

(12d)计算步骤(12c)得到的上包络线和下包络线的平均值m_j-1(t)；

(12e)做减法：h_j(t)＝h_j-1(t)-m_j-1(t)；

(12f)如果满足终止条件，则令c_i(t)＝h_j(t)，然后转入步骤(13)；否则，j＝j+1，转到(12b)；

(13)求剩余信号：r_i(t)＝r_i-1(t)-c_i(t)；

(14)判断r_i(t)的极值点的个数是否仍然超过两个，如果是，则i＝i+1，转到步骤(12)，如果否，则分解结束，r_i(t)即为剩余分量。

3.根据权利要求2所述的高精度的时频分析方法，其特征在于：所述步骤(2)是这样实现的：首先进行MP分解，即搜索组成一个IMF分量的时频原子，然后计算每个时频原子的时频谱，最后将所有时频原子的时频谱相叠加便得到该IMF分量的时频谱，具体包括以下步骤：

(21)初始化，y_j(t)＝c_i(t)，令j＝1；

(22)对y_j(t)进行希尔伯特变换，生成复信号，然后分别求解复信号的瞬时振幅、瞬时频率和瞬时相位；其中u_n为瞬时振幅最大值对应的时间；f_n为该时间对应的瞬时频率；φ_n为该时间对应的瞬时相位；

在一组固定的u_n、f_n和φ_n情况下，通过计算最优化公式得到尺度σ_n；其中D＝{g_r(t)}_r∈Γ为时频原子字典，

是函数R⁽ⁿ⁾ _s和

的内积，

(23)利用所述最优化公式在局域内寻找r_n＝{u_n，σ_n，f_n，φ_n}的最优值，即公式

达到最大时，对应的一组参数便是最优的r_n；搜索范围为[r_n-Δr，r_n+Δr]，其中Δr＝(Δu，Δσ，Δf，Δφ)，Δu为时间偏移量、Δσ为尺度偏移量、Δf为频率偏移量、Δφ为相位偏移量；

(24)在上面两步的基础上，便可以预测到最优的时频原子

的振幅为a_n，即

然后计算剩余信号

令j＝j+1，重复迭代步骤(22)到(24)N次，直到剩余量y_n(t)的最大能量小于原始信号最大能量的10％或更低，得到该IMF分量的N个时频原子，y_n(t)可视为IMF分量c_i(t)的噪音部分；

(25)利用WVD分别计算所述N个时频原子的时频谱，然后将所述N个时频原子的时频谱叠加起来，便得到该IMF分量的时频谱；

(26)对每个IMF分量进行步骤(21)至(25)的处理得到每个IMF分量的时频谱。

Images