技术领域
[0001] 本发明涉及信号处理技术领域，尤其涉及一种信号瞬时频率估计方法。
相关背景技术
[0002] 瞬时频率是表述非平稳信号物理实质的一个特征量，已在雷达、声纳、移动通信、旋转机械等技术领域获得了广泛应用，它可以充分地展现出信号频率随着时间变化的规律，以及能够在时频分布中观测到信号能量集中的程度。因此研究信号瞬时频率的估计方法具有重要的理论意义和实用价值。
[0003] 近些年兴起的信号瞬时频率提取方法主要是将一维的频域和时域信号映射到二维时频域平面上，从而获得信号的时频分布，最终在时频域内完成信号的瞬时频率提取。基于以上理论应运而生了一系列瞬时频率提取方法：短时傅里叶变换(Short‑Time Fourier Transform,STFT)，是一种只具备单一分辨率的处理方法，因此它对突变信号的分析存在局限性，不能敏感地反映信号的突变，只适用于对缓变信号的分析；Wigner‑Ville分布(WVD)算法在计算中不加窗操作，避免了时域分辨率和频域分辨率之间的相互牵制，但在分析多分量信号时WVD会受到交叉项的干扰；小波变换解决了STFT分辨率单一、不能分析非平稳信号的问题，但结果过分依赖于小波基函数，自适应性很差。
[0004] 为提高非平稳信号的瞬时频率提取精度和减少交叉项的干扰，Daubechies Ingrid提出基于小波变换的同步压缩变换，在频率方向提高时频能量分布聚集性且支持信号重构，但在高频段分辨率较低。作为该方法的扩展，Gaurav Thakur提出基于短时傅里叶变换的同步压缩变换，使得在低频段、高频段时频能量分布都具有相同的分辨率。然而，此类同步压缩算法主要存在两方面的不足，一是此类算法仅对瞬时频率恒定的纯谐波信号频率提取精度较高，而在瞬时频率变化剧烈时，同步压缩变换的瞬时频率提取精度较低。因此该算法不适用于非纯谐波信号的瞬时频率提取；二是在实际应用中非平稳轴承振动信号通常包含多个分量，它们不仅在时域内会同时出现，而且瞬时频率还有可能在频域内重叠或各分量的瞬时频率相差较小，此时，利用该类同步压缩算法就不能得到最清晰真实的时频信息。
具体实施方式
[0044] 下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0045] 在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。
[0046] 总体的，如图1所示，本发明实施例公开了一种信号瞬时频率估计方法，包括如下步骤：
[0047] 1)：输入信号记为x(t)，初始化i＝1，xi′(t)＝x(t)；
[0048] 2)：首先将信号xi′(t)按下式变换为X′α(u)：
[0049]
[0050] 其中，
[0051] 其次对X′α(u)进行 搜索得到最优阶次
[0052] 3)：对信号xi′(t)进行下式的变换,此时第i分量的能量绝大部分集中在以该分量真实瞬时频率为中心的一个窄带内，而噪声和其他信号均不会呈现出明显的能量聚集；
[0053]
[0054] 其中，σ为选用的高斯窗函数的标准差，
[0055] 4)：将信号xi′(t)的时频分布按下式进行压缩，将发散的时频能量压缩至真实的瞬时频率处；
[0056]
[0057] 其中，
[0058] 5)：对压缩后的时频能量分布在频率方向进行积分如下式所示，重构信号xi′(t)的时域表示
[0059]
[0060] 其中φk为瞬时频率估计曲线，ds为积分区间；
[0061] 6)：i＝i+1，并对余下的其他分量 重复步骤2)到步骤6)，直至xi′+1(t)幅值低于某一预定阈值。
[0062] 下面结合仿真实验对上述方法进行说明：
[0063] 仿真信号分析：
[0064] 考查如下信号：
[0065]
[0066] 其中， f(t)即为瞬时频率，仿真参数如表1所示。
[0067] 表1仿真参数
[0068]
[0069]
[0070] 仿真信号的时域波形如图2所示，采用本发明所述方法得到信号的时频分布如图3所示。由图3可以看出噪声已经被去除，得到了信号的瞬时频率及2倍和3倍频。为了进一步验证本发明所述方法的有效性，对比分析了本发明所提方法与同步压缩(SST)方法和参数化同步压缩(PSST)方法的瞬时频率估计结果如图4所示。通过图4可知，基于本发明所述方法的瞬时频率估计值在真值附近浮动，而SST和PSST的估计值都小于真值精度较低。
[0071] 为定量的分析本发明所述方法的有效性，利用下式所述的Renyi熵指标和瞬时频率估计误差分析本发明所述方法的时频分布能量聚集性和精度，三种方法的定量对比结果如表2所示。
[0072]
[0073] 其中TF(t,f)表示时频分布，β＝3。Renyi熵越小，表示时频分布能量聚集性越好。
[0074]
[0075] 其中，f(t)为真实瞬时频率， 为瞬时频率估计值，m为采样点数，瞬时频率估计误差越小表示算法的瞬时频率估计精度越高。
[0076] 表2不同算法性能指标对比(多个回车，多一空行)
[0077]
[0078]
[0079] 由表2可知本发明所述方法的Renyi熵和瞬时频率估计误差都小于SST和PSST算法，由此可见本发明所述方法算法不仅提高了信号的时频分布聚集性且提高了瞬时故障频率的估计精度。
[0080] 为了说明本发明所述方法对不同加速度的匀变速信号都具有较强的适用性，在本节中改变转速信号的加速度(调频率)，使调频率在‑20～30范围内以步长2变化得到调频率与瞬时频率估计误差和Renyi熵的关系图如图5，6所示。如图5所示，随着调频率绝对值的增加SST和PSST瞬时频率估计误差逐渐增加，而本发明所述方法不仅对瞬时频率估计误差都小于此两种算法，且随着调频率的增加该误差值变换不大。如图6所示，随着调频率绝对值的增加SST和PSST的Renyi熵逐渐增加，时频聚集性降低。而本发明所述方法的Renyi熵随着调频率的增加变化不大，时频聚集性都高于SST和PSST。
[0081] 为了进一步说明本文所述方法的噪声鲁棒性，改变信号的输入信噪比，使输入信噪比在[‑2 10]dB范围内以步长为1dB变化得到瞬时频率估计误差随输入信噪比的变化曲线如图7所示。如图7可知，随着输入信噪比的增加本发明所述方法不仅使瞬时频率估计误差都小于SST和PSST，而且随着信噪比的增加该值变化不大，说明本发明所述方法对噪声具有较强的鲁棒性。
[0082] 实测信号分析：
[0083] 为验证本发明在实际应用中的有效性，采用QPZZ－Ⅱ旋转机械故障模拟试验台，利用振动信号获得匀减速工况下滚动轴承的瞬时故障频率。试验过程中采样频率为25600Hz，采样时间为4s。试验分析时利用加速度传感器采集滚动轴承的振动信号进行分析，其时域波形如图8所示，瞬时故障频率曲线如图9所示，基于本发明所述方法得到时频分布如图10所示。并利用峰值搜索算法得到本文所述方法，SST和PSST这3种方法的瞬时故障频率估计值，与实际值对比如图11所示。由图11所知，本发明所述方法的瞬时频率估计精度更高，毛刺幅值更小，可见本发明所述方法的抗噪性更强。为定量分析本发明所述方法优势，3种算法的Renyi熵和瞬时故障频率估计误差如表3所示。
[0084] 表3不同算法性能指标对比
[0085]
[0086] 通过表3对比可知，利用本发明所述方法获得的时频能量分布Renyi熵最小，时频能量分布聚集性最高；瞬时故障频率估计误差最小，比SST减少了1.1422Hz，比PSST方法减少了0.8735Hz。由此可见本文所述方法不仅具有较高的时频聚集性、瞬时频率估计精度和抗噪性，也可较好的应用到实际信号瞬时频率提取中。
