CN102650658A - 一种时变非平稳信号时频分析方法 - Google Patents

Abstract

本发明公布了一种时变非平稳信号时频分析方法，首先提出了单频率信号和第i阶局部极值的定义，根据定义，对时变非平稳信号进行单频率信号分解，得到若干个单频率分量，然后通过Hilbert变换获得各个分量的瞬时频率和瞬时幅值。进一步，在幅值-频率-时间三维空间中，把幅值绘制在时频平面上，就得到Hilbert幅度谱。将幅度谱对时间进行积分后，得到Hilbert边际谱。本发明具有完全自适应性；同时不受测不准原理的制约，时间和频率都具有较高的分辨率；分解比较彻底，能得到各自独立的单频率分量即信号的本质特征，各个分量之间不存在交叉项；分解过程中采取加速方法，节约计算时间。

Description

一种时变非平稳信号时频分析方法

技术领域

本发明涉及一种对时变非平稳信号进行时频分析的方法。

背景技术

经典的Fourier变换理论建立了信号从时域到频域的变换桥梁。但是，Fourier变换无法表述频率随时间的变化情况。也就是说，Fourier变换无法表述信号的时频局部特征，仅适合于分析平稳信号而不适合于分析具有时变特性的非平稳信号。

为了弥补Fourier变换这一不足，需要寻求一种能对信号实行时频局部化的新方法来处理信号，便引出了直接在时频平面上表征信号的方法，即时频分析方法。目前应用较多且效果较好的时频分析方法为Hilbert-Huang变换。

1998年，美国国家航空航天局(NASA)的Huang等人提出了Hilbert-Huang变换(HHT)。1999年，Huang又对它进行了一些改进。HHT被认为是近年来对以Fourier变换为基础的线性和稳态谱分析方法的一个重大突破，也是最新发展起来的非线性、非平稳信号的时频分析方法。该方法的基本过程是先将时间信号进行经验模态分解(empirical mode decomposition，EMD)，产生一组具有不同特征时间尺度的本征模函数(intrinsic mode function，IMF)，然后再对每一个IMF分别作Hilbert变换，得到各自的瞬时幅值和瞬时频率。把瞬时幅值表示在时间-频率平面上，就得到了Hilbert谱，该谱能够准确地描述信号的能量随时间和频率的变化规律。该方法计算的瞬时频率有很高的时间分辨率和较高的频率分辨率，然而，这里有三个主要缺点：在低频区会产生不合理的频谱特征；IMF的定义缺乏理论依据，且有的IMF不满足单频率分量条件；不能分离出信号的低能量成分。

发明内容

为了更好地对时变非平稳信号进行时频分析，本发明提供一种新的时变非平稳信号时频分析方法。

本发明的时变非平稳信号进行时频分析的方法，首先将任意时变非平稳信号x(t)，t∈[0，T]分解成若干个单频率分量，包括如下步骤：

步骤一：根据第i阶局部极值的定义，确定信号x(t)，t∈[0，T]所包含的单频率分量的数目n；

步骤二：找出信号x(t)，t∈[0，T]所有的局部极大值点和局部极小值点，用三次样条曲线分别拟合信号的上包络线U(t)和下包络线D(t)，上、下包络线的均值即平均包络线m₁₁(t)为：

$m_{11}\left(t\right)=\frac{U\left(t\right)+D\left(t\right)}{2},$

用局部极值点拟合包络线时，由于信号的起点和终点不一定是局部极值点，因此，需要进行端点处理。本方法采用镜像法对局部极值点进行端点延拓，即在局部极值点序列首尾镜像各增加一个点。

步骤三：将平均包络线m₁₁(t)从原信号x(t)中分离出来，即

h₁₁(t)＝x(t)-m₁₁(t)，

步骤四：判断h₁₁(t)是否满足单频率分量条件，并判断余项x(t)-h₁₁(t)是否满足含有n-1个单频率分量条件，若两个条件同时满足，则h₁₁(t)就是第一阶单频率分量，记为c₁(t)；若有一个条件不满足，则将h₁₁(t)看成x(t)，重复步骤二和步骤三，即h₁₁(t)的平均包络线记为m₁₂(t)，将h₁₁(t)减去m₁₂(t)，得到h₁₂(t)

h₁₂(t)＝h₁₁(t)-m₁₂(t)，

重复以上过程k次，有

h_1k(t)＝h_1(k-1)(t)-m_1k(t)，

综上有

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=h_{11}\left(t\right)-m_{12}\left(t\right)\\ \·\\ \·\\ \·\\ h_{1k}\left(t\right)=h_{1(k-1)}\left(t\right)-m_{1k}\left(t\right)\end{array}\right.,$

迭代终止条件为h_1k(t)满足单频率分量条件，并且余项x(t)-h_1k(t)也满足含有n-1个单频率分量条件，则h_1k(t)就是第一阶单频率分量，记为c₁(t)，它包含了信号的最细尺度或最高频成分。

在本步过程中，为将信号最高频分量快速分离出来，先采取如下加速办法：当k为奇数时，m_1k(t)＝λm_1k(t)，式中λ为加速系数，为大于1的正数，λ取值越大，分离速度越快，但λ取值不宜太大，否则会使上述分离方法不收敛。λ取值应使得在任一时刻t，λm_1k(t)不能超出由其对应的上包络线U(t)和下包络线D(t)确定的区间。为防止λm_1k(t)超出范围致使分离不收敛，另一方面，为加快分离速度，λ取值又较大(比如λ＝5)，因此，可以采取当k大于某一数值(比如k＞10)时开始加速的方法。当k为偶数时，m_1k(t)保持不变。

步骤五：用x(t)减去c₁(t)得到一个去掉高频成分的新信号r₁(t)

r₁(t)＝x(t)-c₁(t)，

将r₁(t)看作x(t)，重复步骤二～步骤四，可以依次得到其它阶单频率分量c₂(t)，c₃(t)，...，c_n(t)和残余项r_n(t)，r_n(t)含有0个单频率分量，即不含单频率分量，代表了信号x(t)的均值或趋势。

由此，x(t)可表示为n个单频率分量和一个残余项的和

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_n\left(t\right),$

步骤六：在单频率信号分解得到了信号x(t)，t∈[0，T]的所有单频率分量以后，分别对每一阶单频率分量应用Hilbert变换，其中第j阶单频率分量c_j(t)的Hilbert变换为

${\hat{c}}_j\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{\∞}\frac{c_j\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ},$

由其构成的解析信号z_j(t)为

$z_j\left(t\right)=c_j\left(t\right)+i{\hat{c}}_j\left(t\right)=A_j\left(t\right)e^{i{\mathrm{\θ}}_j\left(t\right)},$

式中A_j(t)为瞬时幅值，即

$A_j\left(t\right)=\sqrt{c_j^2\left(t\right)+{\hat{c}}_j^2\left(t\right)},$

$\tan \left[{\mathrm{\θ}}_j\left(t\right)\right]=\frac{{\hat{c}}_j\left(t\right)}{c_j\left(t\right)},$

而瞬时频率ω_j(t)为

${\mathrm{\ω}}_j\left(t\right)=\frac{d{\mathrm{\θ}}_j\left(t\right)}{\mathrm{dt}},$

则x(t)可表示为如下形式

$x\left(t\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{c}_j\left(t\right)+r_n\left(t\right)$

$=\mathrm{Re}\left[\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_j\left(t\right)e^{i{\mathrm{\θ}}_j\left(t\right)}\right]+r_n\left(t\right),$

$=\mathrm{Re}\left[\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_j\left(t\right)e^{i{\∫\mathrm{\ω}}_j\left(t\right)\mathrm{dt}}\right]+r_n\left(t\right)$

Re[·]表示取实部。

在振幅-频率-时间三维空间中，把振幅绘制在时频平面上，就得到了Hilbert幅度谱，简称Hilbert谱，它精确地描述了信号的幅值在整个频率段上随频率和时间的变化。由于振幅的平方通常代表能量密度，因此幅度谱也反映了信号能量在空间(或时间)各种尺度上的分布规律。

将幅度谱H(ω，t)对时间t进行积分后，可得到Hilbert边际谱H(ω)

$H\left(\mathrm{\ω}\right)={\∫}_0^{T}H(\mathrm{\ω},t)\mathrm{dt},$

式中，T为信号的采样时间总长度。

把由振幅、频率和时间构成的三维谱图对时间积分，便形成了只有振幅和频率的二维谱图。Hilbert边际谱表征了整个信号中对应每一频率的累积振幅(或能量)的分布情况，它具有一定的概率意义，即在某一频率上存在能量仅意味着具有该频率的波在信号的整个时间跨度内出现的可能性较高。而对于Fourier谱，在某一频率上存在着能量表示的是整个信号里含有该频率的正弦和余弦波。所以，Hilbert边际谱和Fourier谱的意义是不同的。事实上，Hilbert幅度谱可看作是一种加权的非标准化的联合振幅-频率-时间分布，局部振幅就是分配到每个频率-时间单元的权重。从而，Hilbert边际谱中的频率只表明具有该频率的振动有存在的可能性，其确切的发生时间由Hilbert幅度谱给出。

本发明的本质是将信号自适应地由高频到低频进行单频率分量分解，得到若干个单频率分量，通过Hilbert变换获得各个分量的瞬时频率和瞬时幅值，进一步得到幅度谱和边际谱。

本发明不仅可以用于时不变平稳信号，还可以用于时变非平稳信号。本方法可以应用于如下方面：建筑结构、机械振动信号分析，如滤波去噪，模态识别，损伤识别，实时监测等；图像处理；音乐与语言的人工合成；医学成像与诊断；地震勘探数据处理；大型机械的故障诊断等。

单频率分量分解是建立在下列假设基础之上的：

(1)信号至少有两个局部极值点，即一个局部极大值点和一个局部极小值点；

(2)特征时间尺度由连续局部极值点的时间间隔决定；

(3)如果信号没有局部极值点而只有拐点，那么可对其进行一阶或多阶的微分获得局部极值点，最终再通过积分得到相应的分量。

(4)信号需有较高的采样频率，采样频率一般应大于信号最高频率的5倍，即过采样能提高本方法的精度。如信号本身采样频率不够，可以采用插值的方法提高采样频率。

(5)信号的任意两个单频率分量，在任意时刻，两个分量的瞬时频率之比不小于1.1，否则会产生混叠。

本发明具有以下优点：具有完全自适应性；不受测不准原理的制约，时间和频率都具有较高的分辨率；分解比较彻底，能得到各自独立的单频率分量即信号的本质特征，各个分量之间不存在交叉项；分解过程中采取加速方法，可将计算时间缩短为原来的三分之一左右。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书和权利要求书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述，其中：

图1是单频率信号分解的流程图。

图2是本发明的流程图。

图3是考察信号的时域波形图。

图4是考察信号的分解结果图。

图5是考察信号的Hilbert幅度谱图。

图6是考察信号的Hilbert边际谱图。

具体实施方式

以下将参照附图，对本发明的优选实施例进行详细的描述。应当理解，优选实施例仅为了说明本发明，而不是为了限制本发明的保护范围。

为便于阐述本发明，首先提出单频率信号和第i阶局部极值的定义。

1.单频率信号的定义

我们知道，只有对于单频率信号，求解瞬时幅值和瞬时频率才具有实际意义。因此，单频率信号的准确定义至关重要。

以拍振信号为例，如式(1)所示，该信号为调幅-调频信号，其中调幅部分和调频部分均为余弦函数。

y(t)＝2cos(2πf₂t)cos(2πf₁t)，t∈[0，T]，f₁＞10f₂                    (1)

式(1)可变换成式(2)

y(t)＝cos(2π(f₁+f₂)t)+cos(2π(f₁-f₂)t)，t∈[0，T]，f₁＞10f₂          (2)

可以看出，该信号不是单频率信号，而是含有f₁+f₂和f₁-f₂两种频率成分的信号。周期性调幅的本质是信号含有两种不同的频率成分，即信号不是单频率信号。

鉴于此，形如a(t)cos(θ(t))，t∈[0，T]的信号成为单频率信号，需满足如下条件：

①信号的过零点数和局部极值点数必须相等或至多相差1个；

②在信号的任意一点，分别由局部极大值点和局部极小值点确定的上、下包络的均值为0；

③调幅函数a(t)为非周期函数。

第①、②条比较直观，第③条非周期函数比较抽象，不好限定，因此，将调幅函数a(t)成为非周期函数的条件具体化，即需满足如下条件之一：

①调幅函数a(t)不存在局部极值点；

②调幅函数a(t)至多存在一个局部极值点；

③调幅函数a(t)存在若干个局部极值点，但是各个局部极值点的数值几乎相等。

第①条即调幅函数a(t)为单调函数；第②条即调幅函数a(t)至多存在一个局部极大值点或局部极小值点，如果对信号的低频部分不十分关注，可将局部极值点个数加大；第③条即调幅函数a(t)为常函数，理论上不存在局部极值点，但是由于采样、计算等误差的影响，出现局部极值点，但是数值几乎相等。

对于任意时变非平稳信号x(t)，可以表示为如式(3)所示的若干个单频率分量之和的形式：

$x\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_i\left(t\right)\cos \left({\mathrm{\θ}}_i\left(t\right)\right),t\∈[0,T]---\left(3\right)$

式中，a_i(t)为第i个单频率分量的幅值调制函数，为非周期函数；θ_i(t)为第i个单频率分量的相位调制函数。

2.第i阶局部极值的定义

时频分析的目的是得到瞬时频率和瞬时幅值。信号的本质就是波动，其局部极值无疑能反映信号的特性，如频率和幅值。对任意信号x(t)，t∈[0，T]，局部极大值集为：

$S_{0_{\max }}=\left\{s_{0_{{\max ,}_{t_k}}}\right|s_{0_{{\max ,}_{t_k}}}=x\left(t_k\right)\&cap;x^{\&prime;}\left(t_k\right)=0\&cap;x^{\&prime;\&prime;}\left(t_k\right)<0\},t_k\&Element;[0,T]---\left(4\right)$

局部极小值集为：

$S_{0_{\min }}=\left\{s_{0_{{\min ,}_{t_k}}}\right|s_{0_{{\min ,}_{t_k}}}=x\left(t_k\right)\&cap;x^{\&prime;}\left(t_k\right)=0\&cap;x^{\&prime;\&prime;}\left(t_k\right)>0\},t_k\&Element;[0,T]---\left(5\right)$

定义：

第0阶局部极大极大值集为：

$S_{0_{\max \max }}=S_{0_{\max }}---\left(6\right)$

第0阶局部极大极大值集为：

$S_{0_{\max \min }}=S_{0_{\max }}---\left(7\right)$

第0阶局部极小极大值集为：

$S_{0_{\min \max }}=S_{0_{\min }}---\left(8\right)$

第0阶局部极小极大值集为：

$S_{0_{\min \min }}=S_{0_{\min }}---\left(9\right)$

第i(i＝1，2，3，…)阶局部极大极大值集为：

$S_{i_{\max \max }}=\left\{s_{i_{{\max \max ,}_{t_k}}}\right|s_{i_{{\max \max ,}_{t_k}}}=s_{{i-1}_{{\max \max ,}_{t_k}}}\&cap;s_{{i-1}_{{\max \max }_{t_k}}}>s_{{i-1}_{{\max \max }_{t_{k-1}}}}---\left(10\right)$

$\&cap;s_{{i-1}_{{\max \max }_{t_k}}}>s_{{i-1}_{{\max \max }_{t_{k+1}}}}\},t_{k-1},t_k,t_{k+1}\&Element;[0,T]$

第i阶局部极大极小值集为：

$S_{i_{\max \min }}=\left\{s_{i_{\max \min {,}_{t_k}}}\right|s_{i_{{\max \min ,}_{t_k}}}=s_{{i-1}_{{\max \max ,}_{t_k}}}\&cap;s_{{i-1}_{{\max \max }_{t_k}}}<s_{{i-1}_{\max {\max }_{t_{k-1}}}}---\left(11\right)$

$\&cap;s_{{i-1}_{{\max \max }_{t_k}}}<s_{{i-1}_{{\max \max }_{t_{k+1}}}}\},t_{k-1},t_k,t_{k+1}\&Element;[0,T]$

第i阶局部极小极大值集为：

$S_{i_{\min \max }}=\left\{s_{i_{\min \max {,}_{t_k}}}\right|s_{i_{{\min \max ,}_{t_k}}}=s_{{i-1}_{\min \min {,}_{t_k}}}\&cap;s_{{i-1}_{{\min \min }_{t_k}}}>s_{{i-1}_{{\min \min }_{t_{k-1}}}}---\left(12\right)$

$\&cap;s_{{i-1}_{{\max \min }_{t_k}}}>s_{{i-1}_{{\max \min }_{t_{k+1}}}}\},t_{k-1},t_k{,t}_{k+1}\&Element;[0,T]$

第i阶局部极小极小值集为：

$S_{i_{\min \min }}=\left\{s_{i_{{\min \min ,}_{t_k}}}\right|s_{i_{\min \min {,}_{t_k}}}=s_{{i-1}_{{\min \min ,}_{t_k}}}\&cap;s_{{i-1}_{{\min \min }_{t_k}}}<s_{{i-1}_{{\min \min }_{t_{k-1}}}}---\left(13\right)$

若信号x(t)，t∈[0，T]的第i阶局部极大极大值集

(或第i阶局部极小极小值集

满足如下条件之一，则该信号含有i个单频率分量。

①

或

为空集但

或

均至少有2个元素；

②

或

只含有一个元素；

③

或

含有若干个元素，但是各个元素的数值几乎相等。

第①条即

或

为单调序列；第②条即或

至多存在一个局部极大值点或局部极小值点，如果不十分关注信号的低频部分，可将元素个数加大；第③条即或

为常数序列，理论上

或

为空集，但是由于采样、计算等误差的影响，使得

或

为非空集，但是各个元素数值几乎相等。

根据上述定义，本发明的一种时变非平稳信号时频分析方法，包括以下步骤：

步骤一：输入原始信号至处理机，根据第i阶局部极值的定义，确定信号x(t)，t∈[0，T]所包含的单频率分量的数目n；

步骤二：找出信号x(t)，t∈[0，T]所有的局部极大值点和局部极小值点，用三次样条曲线分别拟合信号的上包络线U(t)和下包络线D(t)，上、下包络线的均值即平均包络线m₁₁(t)为：

$m_{11}\left(t\right)=\frac{U\left(t\right)+D\left(t\right)}{2}---\left(14\right)$

用局部极值点拟合包络线时，由于信号的起点和终点不一定是局部极值点，因此，需要进行端点处理。本方法采用镜像法对局部极值点进行端点延拓，即在局部极值点序列首尾镜像各增加一个点；

步骤三：将平均包络线m₁₁(t)从原信号x(t)中分离出来，即

h₁₁(t)＝x(t)-m₁₁(t)                                       (15)

步骤四：判断h₁₁(t)是否满足单频率分量条件，并判断余项x(t)-h₁₁(t)是否满足含有n-1个单频率分量条件，若两个条件同时满足，则h₁₁(t)就是第一阶单频率分量，记为c₁(t)；若有一个条件不满足，则将h₁₁(t)看成x(t)，重复步骤二和步骤三，即h₁₁(t)的平均包络线记为m₁₂(t)，将h₁₁(t)减去m₁₂(t)，得到h₁₂(t)

h₁₂(t)＝h₁₁(t)-m₁₂(t)                                     (16)

重复以上过程k次，有

h_1k(t)＝h_1(k-1)(t)-m_1k(t)                                 (17)

综上有

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=h_{11}\left(t\right)-m_{12}\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ h_{1k}\left(t\right)=h_{1(k-1)}\left(t\right)-m_{1k}\left(t\right)\end{array}\right.---\left(18\right)$

迭代终止条件为h_1k(t)满足单频率分量条件，并且余项x(t)-h_1k(t)也满足含有n-1个单频率分量条件，则h_1k(t)就是第一阶单频率分量，记为c₁(t)，它包含了信号的最细尺度或最高频成分。

在本步过程中，为将信号最高频分量快速分离出来，先采取如下加速办法：当k为奇数时，m_1k(t)＝λm_1k(t)，式中λ为加速系数，为大于1的正数，λ取值越大，分离速度越快，但λ取值不宜太大，否则会使上述分离方法不收敛。λ取值应使得在任一时刻t，λm_1k(t)不能超出由其对应的上包络线U(t)和下包络线D(t)确定的区间。为防止λm_1k(t)超出范围致使分离不收敛，另一方面，为加快分离速度，λ取值又较大(比如λ＝5)，因此，可以采取当k大于某一数值(比如k＞10)时开始加速的方法。当k为偶数时，m_1k(t)保持不变。

步骤五：用x(t)减去c₁(t)得到一个去掉高频成分的新信号r₁(t)

r₁(t)＝x(t)一c₁(t)                                          (19)

将r₁(t)看作x(t)，重复步骤二～步骤四，可以依次得到其它阶单频率分量c₂(t)，c₃(t)，...，c_n(t)和残余项r_n(t)，r_n(t)含有0个单频率分量，即不含单频率分量，代表了信号x(t)的均值或趋势。

由此，x(t)可表示为n个单频率分量和一个残余项的和

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right)---\left(20\right)$

本分解方法是基于单频率信号的定义进行的，因此，将本分解方法命名为单频率信号分解。上述步骤一～步骤五可以视为单频率信号分解的内容。

步骤六：在单频率信号分解得到了信号x(t)，t∈[0，T]的所有单频率分量以后，分别对每一阶单频率分量应用Hilbert变换，其中第j阶单频率分量c_j(t)的Hilbert变换为

${\hat{c}}_j\left(t\right)=\frac{1}{\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{\&infin;}\frac{c_j\left(\mathrm{\&tau;}\right)}{t-\mathrm{\&tau;}}\mathrm{d\&tau;}---\left(21\right)$

由其构成的解析信号z_j(t)为

$z_j\left(t\right)=c_j\left(t\right)+i{\hat{c}}_j\left(t\right)=A_j\left(t\right)e^{i{\mathrm{\&theta;}}_j\left(t\right)}---\left(22\right)$

式中A_j(t)为瞬时幅值，即

$A_j\left(t\right)=\sqrt{c_j^2\left(t\right)+{\hat{c}}_j^2\left(t\right)}---\left(23\right)$

$\tan \left[{\mathrm{\&theta;}}_j\left(t\right)\right]=\frac{{\hat{c}}_j\left(t\right)}{c_j\left(t\right)}---\left(24\right)$

而瞬时频率ω_j(t)为

${\mathrm{\&omega;}}_j\left(t\right)=\frac{d{\mathrm{\&theta;}}_j\left(t\right)}{\mathrm{dt}}---\left(25\right)$

则x(t)可表示为如下形式

$x\left(t\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_j\left(t\right)+r_n\left(t\right)$

$=\mathrm{Re}\left[\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_j\left(t\right)e^{i{\mathrm{\&theta;}}_j\left(t\right)}\right]+r_n\left(t\right)---\left(26\right)$

$=\mathrm{Re}\left[\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_j\left(t\right)e^{i{\&Integral;\mathrm{\&omega;}}_j\left(t\right)\mathrm{dt}}\right]+r_n\left(t\right)$

Re[·]表示取实部。

在振幅-频率-时间三维空间中，把振幅绘制在时频平面上，就得到了Hilbert幅度谱，简称Hilbert谱，它精确地描述了信号的幅值在整个频率段上随频率和时间的变化。由于振幅的平方通常代表能量密度，因此幅度谱也反映了信号能量在空间(或时间)各种尺度上的分布规律。

将幅度谱H(ω，t)对时间t进行积分后，可得到Hilbert边际谱H(ω)

$H\left(\mathrm{\&omega;}\right)={\&Integral;}_0^{T}H(\mathrm{\&omega;},t)\mathrm{dt}---\left(27\right)$

式中，T为信号的采样时间总长度。

实施例

考察如式(28)所示的信号：

x(t)＝2cos(20πt)cos[100πt+2cos(6πt)]+2e^0.2tsin(10πt+2πt²)               (28)

采样频率取1000Hz，采样时间取0～2s，其时域波形如图3所示。

式(28)可变为

x(t)＝cos[120πt+2cos(6πt)]+cos[80πt+2cos(6πt)]+2e^0.2tsin(10πt+2πt²)    (29)

采用本发明方法对其进行单频率信号分解，分解得到c1、c2、c 3分量和趋势项r(t)，如图4所示。由图4可以看出，三个分量分别对应仿真信号x(t)的三个分量，因此，本发明方法是根据信号本身而进行的自适应分解，得到的每-个分量都具有一定的物理意义，反映了信号的内在本质。在振幅-频率-时间三维空间中，把振幅绘制在时频平面上，就得到了Hilbert幅度谱，如图5所示，图中色条颜色表示幅值大小。将Hilbert幅度谱对时间进行积分后，可得到Hilbert边际谱，如图6所示。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种时变非平稳信号时频分析方法，其特征在于：包括以下步骤：

步骤一：输入原始信号，根据第i阶局部极值的定义，确定信号x(t)，t∈[0，T]所包含的单频率分量的数目n；

步骤二：找出信号x(t)，t∈[0，T]所有的局部极大值点和局部极小值点，用三次样条曲线分别拟合信号的上包络线U(t)和下包络线D(t)，上、下包络线的均值即平均包络线m₁₁(t)为：

$m_{11}\left(t\right)=\frac{U\left(t\right)+D\left(t\right)}{2},$

步骤三：将平均包络线m₁₁(t)从原信号x(t)中分离出来，即

h₁₁(t)＝x(t)-m₁₁(t)，

步骤四：判断h₁₁(t)是否满足单频率分量条件，并判断余项x(t)-h₁₁(t)是否满足含有n-1个单频率分量条件，若两个条件同时满足，则h₁₁(t)就是第一阶单频率分量，记为c₁(t)；若有一个条件不满足，则将h₁₁(t)看成x(t)，重复步骤二和步骤三，即h₁₁(t)的平均包络线记为m₁₂(t)，将h₁₁(t)减去m₁₂(t)，得到h₁₂(t)

h₁₂(t)＝h₁₁(t)-m₁₂(t)，

重复以上过程k次，有

h_1k(t)＝h_1(k-1)(t)-m_1k(t)，

综上有

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=h_{11}\left(t\right)-m_{12}\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ h_{1k}\left(t\right)=h_{1(k-1)}\left(t\right)-m_{1k}\left(t\right)\end{array}\right.,$

迭代终止条件为h_1k(t)满足单频率分量条件，并且余项x(t)-h_1k(t)也满足含有n-1个单频率分量条件，则h_1k(t)就是第一阶单频率分量，记为c₁(t)，它包含了信号的最细尺度或最高频成分；

步骤五：用x(t)减去c₁(t)得到一个去掉高频成分的新信号r₁(t)

r₁(t)＝x(t)-c₁(t)，

将r₁(t)看作x(t)，重复步骤二～步骤四，可以依次得到其它阶单频率分量c₂(t)，c₃(t)，...，c_n(t)和残余项r_n(t)，r_n(t)含有0个单频率分量，即不含单频率分量，代表了信号x(t)的均值或趋势；

由此，x(t)可表示为n个单频率分量和一个残余项的和

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right);$

步骤六：在单频率信号分解得到了信号x(t)，t∈[0，T]的所有单频率分量以后，分别对每一阶单频率分量应用Hilbert变换，其中第j阶单频率分量c_j(t)的Hilbert变换为

${\hat{c}}_j\left(t\right)=\frac{1}{\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{\&infin;}\frac{c_j\left(\mathrm{\&tau;}\right)}{t-\mathrm{\&tau;}}\mathrm{d\&tau;},$

由其构成的解析信号z_j(t)为

$z_j\left(t\right)=c_j\left(t\right)+i{\hat{c}}_j\left(t\right)=A_j\left(t\right)e^{i{\mathrm{\&theta;}}_j\left(t\right)},$

式中A_j(t)为瞬时幅值，即

$A_j\left(t\right)=\sqrt{c_j^2\left(t\right)+{\hat{c}}_j^2\left(t\right)},$

$\tan \left[{\mathrm{\&theta;}}_j\left(t\right)\right]=\frac{{\hat{c}}_j\left(t\right)}{c_j\left(t\right)},$

而瞬时频率ω_j(t)为

${\mathrm{\&omega;}}_j\left(t\right)=\frac{d{\mathrm{\&theta;}}_j\left(t\right)}{\mathrm{dt}},$

则x(t)可表示为如下形式

$x\left(t\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_j\left(t\right)+r_n\left(t\right)$

$=\mathrm{Re}\left[\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_j\left(t\right)e^{i{\mathrm{\&theta;}}_j\left(t\right)}\right]+r_n\left(t\right),$

$=\mathrm{Re}\left[\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_j\left(t\right)e^{i{\&Integral;\mathrm{\&omega;}}_j\left(t\right)\mathrm{dt}}\right]+r_n\left(t\right)$

Re[·]表示取实部；

在振幅-频率-时间三维空间中，把振幅绘制在时频平面上，就得到了Hilbert幅度谱，简称Hilbert谱；

将幅度谱H(ω，t)对时间t进行积分后，可得到Hilbert边际谱H(ω)

$H\left(\mathrm{\&omega;}\right)={\&Integral;}_0^{T}H(\mathrm{\&omega;},t)\mathrm{dt},$

式中，T为信号的采样时间总长度。

2.根据权利要求1所述的一种时变非平稳信号时频分析方法，其特征在于：在步骤二中，采用镜像法对局部极值点进行端点延拓，即在局部极值点序列首尾镜像各增加一个点。

3.根据权利要求1所述的一种时变非平稳信号时频分析方法，其特征在于：在步骤四中，为将信号最高频分量快速分离出来，先采取如下加速办法：当k为奇数时，m_1k(t)＝λm_1k(t)，式中λ为加速系数，为大于1的正数，λ取值应使得在任一时刻t，λm_1k(t)不能超出由其对应的上包络线U(t)和下包络线D(t)确定的区间。

Images