CN103020479A - 基于非线性调频小波变换的信号瞬时频率估计方法 - Google Patents

Abstract

本发明公开了一种估计非平稳信号瞬时频率的方法，用于信号处理领域。本发明利用非线性调频小波变换技术，考虑到非线性调频小波变换中多项式核特征参数对时频分析性能的影响，通过选择合适的核特征参数，实现对非平稳信号瞬时频率的精确估计。本发明的优点在于：1.由于采用非线性调频小波变换技术，从而达到精确估计非平稳信号瞬时频率的目的，其关键技术是非线性调频小波变换的非线性多项式核特征参数选择。2.利用基于非线性调频小波变换的瞬时频率估计技术，可以避免常规技术上计算繁杂等缺点，具有实现简单，精度高的优点。3.在算法上简单易行，可适用信号处理等技术领域。

Description

基于非线性调频小波变换的信号瞬时频率估计方法

技术领域

本发明涉及信号处理领域，尤其是一种估计非平稳信号瞬时频率的方法，具体是一种基于非线性调频小波变换的瞬时频率估计方法。

背景技术

在实际工程应用场合，非平稳信号是普遍存在的，如语音、生物医学、雷达、声纳和地震监测信号及各种机械振动等，其共同特征是存在着变化的频率，即其统计特性(相关函数、功率谱等)是时变函数。传统的傅立叶变换(FT)针对的是周期性平稳信号，它依赖于信号的全局信息，并不能反映信号的局部特征，因此，对分析非平稳信号不具有效性。而瞬时频率IF(Instantaneous Frequency)在描绘非平稳信号的特征时具有它独特的瞬时有效性。

在众多瞬时频率估计（IF Estimation）方法中，使用基于时频分布（TFD）的方法得到的结果更为可靠并且具有更好的抗噪性能，因而得到了广泛的关注与研究。基于时频分布的方法利用了信号的时频分布在瞬时频率或其附近能量集中的特性。经对现有技术文献的检索发现，三种常用来做时频分布的方法分别是：短时傅里叶变换（STFT）、连续小波变换（CWT）和Wigner-Ville分布（WVD）。短时傅里叶变换和连续小波变换是一种线性时频分析方法，在时频平面有着静态的分辨率。由于Heisenberg测不准原理的限制，这两种方法不能同时取得很好的时域分辨率和频域分辨率，只能实现有限精度的瞬时频率估计，尤其是当瞬时频率随时间做非线性变化时。Wigner-Ville分布(WVD)是一种双线性变换，有着良好的数学特性，可以精确地反映线性调频信号的频率信息。但遗憾的是，这类分布对频率随时间呈非线性变化的信号或包含多个分量的信号，不可避免地存在交叉项，使得很难识别自身项信号。很多学者针对如何抑制干扰项的问题做了大量工作：Stankovic和他的同事们提出了一种自适应窗长度的WVD方法和变化的数据驱动的WVD方法；Boashash等人提出了高阶Wigner-Ville分布(PWVD)；还有通过加窗处理得到的伪WVD和平滑伪WVD。但这些方法要么缺乏普遍适用性，要么是以牺牲时频聚集性为代价的。

除了以上提到的时频分析方法以外，线性调频小波变换（Chirplet）是一种特别针对线性调频信号的时频分析方法。针对线性调频小波变换中的多参数导致的计算过程繁杂问题，O’Neil和Flandrin等人提出了自适应Chirplet分解方法，通过求解最佳的线性调频小波参数以便对信号进行匹配。但当研究信号的频率随时间呈强烈非线性变化时，线性调频小波时频分布的聚集性很差，不能保证其分析的精确性。

本发明提供一种新的参数化时频分析方法-非线性调频小波变换（Nonlinear ChirpletTransform）。基于非线性调频小波变换发展出了一种针对严重非线性频率族信号的瞬时频率估计方法。非线性调频小波变换用新的多项式非线性瞬时频率律核函数取代了线性瞬时频率律的Chirplet核函数。Weierstrass逼近定理表明闭区间上的连续函数可用多项式级数一致逼近。因此，对瞬时频率是时间任意连续函数的信号而言，非线性调频小波变换是一种有着更高精度的分析方法。至今为止，通过大量的文献检索，还没有发现采用非线性调频小波变换技术估计信号瞬时频率的任何报道。

发明内容

本发明的目的在于克服现有技术中存在的分析瞬时频率是时间的非线性函数信号的不足和缺陷，提供一种基于非线性调频小波变换的瞬时频率估计方法,通过多项式非线性核函数取代线性调频小波变换中的线性调频核函数，非线性调频小波变换的时频分布结果有良好的时频聚集性，且容易识别出信号的时频特征。

根据本发明的一个方面，提供一种基于非线性调频小波变换的信号瞬时频率估计方法，首先利用短时傅里叶变换得到初始化的时频分布，沿频率轴方向提取时频分布的局部极大值曲线，以得到信号瞬时频率的粗略估计；然后利用最小二乘法拟合逼近该粗略估计的瞬时频率，将拟合多项式系数就作为多项式核特征参数对信号做非线性调频小波变换；最后，在信号的非线性调频小波变换结果上，再次提取时频分布的局部极大值曲线，将该再次提取时频分布的局部极大值曲线随时间的变化规律作为该信号的瞬时频率估计。

优选地，具体包括如下步骤：

第1步：设终止条件阈值δ；令多项式核函数的特征参数(α₁,…,α_n)＝0；设定高斯窗口长度；

第k步，k＞1：具体包括如下子步骤：

1）作特征参数为(α₁,…,α_n)的非线性调频小波变换；

2）提取时频分布的极大值曲线得到信号的瞬时频率估计

3）采用最小二乘法用n阶多项式逼近信号的瞬时频率估计

记估计得到的系

数为

4）计算终止判据ξ_(s)；

5）如果ξ_(s)＞δ，则令

k＝k+1，转到步骤1)继续执行；否

则转到步骤6）；

6）取关于时间的n阶多项式

作为信号的瞬时频率估计；

7）退出。

更为具体地，本发明是通过以下技术方案实现的，利用新的参数化时频分析方法-非线性调频小波变换，考虑到非线性调频小波变换中多项式核特征参数对时频分析性能的影响，通过选择合适的核特征参数，实现对非平稳信号瞬时频率的精确估计。本发明提出的非平稳信号瞬时频率估计方法具有实现简单，精度高的优点。

以下对本发明方法作进一步的说明，包括非线性调频小波变换的基本原理和基于非线性调频小波变换的瞬时频率估计方法两个部分的介绍。

1.非线性调频小波变换

对于信号s(t)∈L²(R)，其非线性调频小波变换定义为

${\mathrm{NCT}}_s(t_0,\mathrm{\ω},{\mathrm{\α}}_1,...,{\mathrm{\α}}_n;\mathrm{\σ})={\∫}_{-\∞}^{+\∞}z\left(t\right){\mathrm{\Φ}}_{{\mathrm{\α}}_1,...,{\mathrm{\α}}_n}^R\left(t\right){\mathrm{\Φ}}_{{\mathrm{\α}}_1,...,{\mathrm{\α}}_n}^M(t,t_0)w_{\left(\mathrm{\σ}\right)}(t-t_0)\exp (-\mathrm{j\ωt})\mathrm{dt}$ 其中，

$\left\{\begin{array}{c}\stackrel{\‾}{z}\left(t\right)=z\left(t\right){\mathrm{\Φ}}_{\mathrm{\α}}^R\left(t\right){\mathrm{\Φ}}_{\mathrm{\α}}^M(t,t_0)\\ {\mathrm{\Φ}}_{{\mathrm{\α}}_1,...,{\mathrm{\α}}_n}^R\left(t\right)=\exp (-j\underset{k=2}{\overset{n+1}{\mathrm{\Σ}}}\frac{1}{k}{\mathrm{\α}}_{k-1}{t}^k)\\ {\mathrm{\Φ}}_{{\mathrm{\α}}_1,...,{\mathrm{\α}}_n}^M(t,t_0)=\exp \left(j\underset{k=2}{\overset{n+1}{\mathrm{\Σ}}}{\mathrm{\α}}_{k-1}{t}_0^{(k-1)}t\right)\end{array}\right.$

z(t)是信号s(t)的解析信号，由Hilbert变换得到，z(t)＝s(t)+jH[s(t)]；参数t₀∈R代表时间；w∈L²(R)代表一个非负的对称的归一化的实窗函数，通常取高斯函数，表达式为 $w_{\left(\mathrm{\σ}\right)}\left(t\right)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp (-\frac{1}{2}{(t/\mathrm{\σ})}^2),$

和

分别是多项式频率旋转算子和频率平移算子，(α₁,…,α_n)是多项式核函数的特征参数。

非线性调频小波变换的运算原理如图1所示，其中f_s(t)是信号的瞬时频率。具体地，首先将信号在时频面内旋转，即将信号的瞬时频率f_s(t)减去非线性调频核函数的瞬时频率

然后以

为增量平移频率，再做窗函数为w_(σ)的短时傅里叶变换。显然，在特殊时刻t₀处，非线性调频小波变换的频率分辨率由两部分决定：一是时间段[t₀-σ/2,t₀+σ/2]内

的频率范围，记为Δf_s(t₀;σ)，二是高斯窗的频宽1/σ；即t₀时刻的频率分辨率为f_s(t₀;σ)+1/σ。理想情况下，如果非线性调频核函数的瞬时频率与信号的频率曲线正好相吻合，那么为一常数ω₀，Δf_s(t₀;σ)在该时间段内的频率范围为零，因此非线性调频小波变换在该时刻的频率分辨率达到最小，为1/σ。

2.基于非线性调频小波变换的瞬时频率估计方法

本发明中，基于非线性调频小波变换的瞬时频率估计方法的基本思想是：首先利用短时傅里叶变换得到初始化的时频分布，沿频率轴方向提取时频分布的局部极大值曲线以得到信号瞬时频率的粗估计；然后利用最小二乘法拟合逼近该粗略估计的瞬时频率，将拟合多项式系数就作为多项式核特征参数对信号做非线性调频小波变换；最后，信号频谱分量上能量峰值的时间分布就作为该信号的瞬时频率估计。由前面的介绍可以知道，对非线性调频小波变换来说，当其非线性调频核函数的瞬时频率与信号的瞬时频率越一致，其时频分布的聚集性会越好。因此，虽然用短时傅立叶变换结果粗略估计得到的瞬时频率与真实瞬时频率可能并不十分一致，但其后的非线性调频小波变换结果将会更接近真实的瞬时频率。因此可用非线性调频小波变换得到的瞬时频率更新多项式核，然后用新的多项式核再做变换。这个过程可以一直重复进行，直到估计得到的瞬时频率没有明显的变化。该过程的终否，可简单地用连续两次估计得到的瞬时频率即IF_(i+1)(t)和IF_(i)(t)之间的差异作为判断条件，终止判据如下

${\mathrm{\&xi;}}_{\left(s\right)}=\mathrm{mean}(\&Integral;\frac{|{\mathrm{IF}}_{(i+1)}\left(t\right)-{\mathrm{IF}}_{\left(i\right)}\left(t\right)|}{\left|{\mathrm{IF}}_{\left(i\right)}\left(t\right)\right|}\mathrm{dt})<\mathrm{\&delta;}$

其中，δ为预先设定的阈值。

基于非线性调频小波变换的瞬时频率估计技术可以归纳为如下步骤：

第1步：

设定阈值δ；令(α₁,…,α_n)＝0；设定高斯窗口长度；

第k步：

1)作特征参数为(α₁,…,α_n)的非线性调频小波变换；

2)提取时频分布的极大值曲线得到信号的瞬时频率估计

3)采用最小二乘法用n阶多项式逼近信号的瞬时频率估计

记估计得到的系

数为

4)计算终止判据ξ_(s)；

5)如果ξ_(s)＞δ，

k＝k+1，转到步骤1)；

否则转到步骤6)；

6)取关于时间的n阶多项式

作为信号的瞬时频率估计，；

7)退出

与现有技术相比，本发明的优点在于：

1.由于采用非线性调频小波变换技术，从而达到精确估计非平稳信号瞬时频率的目的，其关键技术是非线性调频小波变换的非线性多项式核特征参数选择。

2.利用基于非线性调频小波变换的瞬时频率估计技术，可以避免常规技术上计算繁杂等缺点，具有实现简单，精度高的优点。

3.在算法上简单易行，可适用信号处理等技术领域。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为非线性调频小波变换的原理图；

图2为仿真Jeffcott转子的振动信号；

图3为第一次非线性调频小波变换后的时频分布和峰值曲线；

图4为第三次非线性调频小波变换后的时频分布和峰值曲线；

图5为转子实验台；

图6为实测转子的振动信号；

图7为第一次非线性调频小波变换后的时频分布和峰值曲线；

图8为第三次非线性调频小波变换后的时频分布和峰值曲线。

图中：1为电动机，2为转子，3为加速度传感器。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

旋转机械经常经历起动、停机和增减负荷等过渡过程，这些过程中的非平稳响应信号蕴含着丰富的机组状态信息。因此，对机组过渡过程中出现的现场测试信号进行分析处理，提取信号在时频域内的特征，对于旋转机械及检测诊断而言具有重要意义。因此，本发明中用非线性调频小波变换来估计起停机过程中转子的瞬态转速。

为更好地理解本发明的技术方案，以下结合附图及具体的实施例作进一步描述。

实施例1

图2为仿真Jeffcott转子启动过程中的振动响应信号，转速呈指数规律升速，采样频率为500Hz。采用非线性调频小波变换算法时，高斯窗长度为512，多项式核函数的阶数设为4阶，终止条件为

${\mathrm{\&xi;}}_{\left(s\right)}=\mathrm{mean}(\&Integral;\frac{|{\mathrm{IF}}_{(i+1)}\left(t\right)-{\mathrm{IF}}_{\left(i\right)}\left(t\right)|}{\left|{\mathrm{IF}}_{\left(i\right)}\left(t\right)\right|}\mathrm{dt})<\mathrm{\&delta;}$

其中阈值δ设为0.1%。当判定条件达到阈值之前，共进行了3次非线性调频小波变换运算；图3和图4分别表示第一次和第三次变换的时频分布结果、提取的时频分布极大值曲线和估计的转速曲线。显然，图4的时频聚集性比图3中短时傅里叶变换的结果要好很多（图中能量高的区域是因为转子在启动过程中达到一阶临界转速引起的共振）。启动过程中转子的瞬态转速可用多项式函数的估计系数表示如下：

IS(t)≈0.469+2.5448t+1.162t²-0.0808t³+0.0161t⁴(Hz)

实施例2

图5为实验装置图，振动信号由加速度传感器测得。图6为实测的振动信号，采样频率为100Hz。采用非线性调频小波变换算法时，高斯窗长度为512，多项式核函数的阶数设为7阶，阈值δ设为0.1%。当判定条件达到阈值之前，共进行了3次非线性调频小波变换运算；图7和图8分别表示第一次和第三次变换的时频分布结果、提取的时频分布极大值曲线和估计的转速曲线。显然，图8的时频聚集性比图7中短时傅里叶变换的结果要好很多。停机过程中转子的瞬态转速可用多项式函数的估计系数表示如下：

IS(t)≈30.6692-0.1367t+2.0530t²-1.6297t³

+0.5266t⁴-0.0905t⁵+0.0090t⁶-0.0005t⁷(Hz)

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (2)

1.一种基于非线性调频小波变换的信号瞬时频率估计方法，其特征在于，首先利用短时傅里叶变换得到初始化的时频分布，沿频率轴方向提取时频分布的局部极大值曲线，以得到信号瞬时频率的粗略估计；然后利用最小二乘法拟合逼近该粗略估计的瞬时频率，将拟合多项式系数就作为多项式核特征参数对信号做非线性调频小波变换；最后，在信号的非线性调频小波变换结果上，再次提取时频分布的局部极大值曲线，将该再次提取时频分布的局部极大值曲线随时间的变化规律作为该信号的瞬时频率估计。

2.根据权利要求1所述的基于非线性调频小波变换的信号瞬时频率估计方法，其特征在于，具体包括如下步骤：

第1步：设定终止条件阈值δ；令多项式核函数的特征参数(α₁,…,α_n)＝0；设定高斯窗口长度；

第k步，k＞1：具体包括如下子步骤：

1)作特征参数为(α₁,…,α_n)的非线性调频小波变换；

2)提取时频分布的极大值曲线得到信号的瞬时频率估计

3)采用最小二乘法用n阶多项式逼近信号的瞬时频率估计

并记估计得到的系数为

4)计算终止判据ξ_(s)；

5)如果ξ_(s)＞δ，则令

k＝k+1，转到步骤1)继续执行；否则转到步骤6）；

6)取关于时间的n阶多项式

作为信号的瞬时频率估计；

7)退出。

Images