频率估计的时频分析方法研究

郭 振，郭建涛

（1.信阳师范学院 计算机与信息技术学院，河南 信阳 464000；2.信阳师范学院 物理电子工程学院，河南 信阳 464000）

信号处理与分析中最重要也是最基本的两个变量是时间和频率。随着现代电子技术的发展使其信号形式越来越复杂，信号的非线性、非平稳时变特性已经成为信号的主要特征，单纯的时域或者频域处理已经不能满足需要。采用时频分析可以在时频二维平面给出信号的时变频谱，是解决上述问题的一种有效方法。

现代通信，特别是跳频扩频通信系统中，通过频率和相位调制来展宽信号的频谱。在雷达、声纳以及电子对抗中的信号探测系统中，系统和目标之间的相对移动，传播媒质的扰动都会使频率发生改变。信号频率的时变已经成为其非平稳性的主要表现之一，因此在很多实际信号的处理中，估计一个非平稳过程的瞬时频率是一项重要的工作。利用时频分析研究瞬时频率估计及各种方法的性能对国民经济的发展和国防建设有重要的意义。笔者阐述了瞬时频率时频分析理论研究和应用发展概况，给出了几种典型的估计算法以及相应的性能评价。

1 瞬时频率

瞬时频率这个概念起源于通信中的调频研究。Armstrong[1]发现进行频率调制可以有效地压制噪声，从而引起人们对调频信号中“频率”这个物理量的兴趣，促使人们去深入研究频率调制和瞬时频率的概念及其数学描述。

1）解析信号的相位差分

Gabor[2]提出了从实信号产生相应复信号的方法，这对瞬时频率概念的进一步发展起了重要作用。实信号s（t）=a（t）cos（φ（t））的瞬时频率定义式为：

其中，z（t）为实信号的解析信号，φ（t）为解析信号相位。

2）基于时频分布矩（特定时间的频率平均）：许多TFD例如WVD的一阶矩给出瞬时频率值，另外一些TFD例如STFT的一阶矩给出瞬时频率的近似值。定义式为：

应用时频分布WVD的时间一阶条件矩计算得到的瞬时频率和应用解析信号有限差分法得到的瞬时频率结果完全相同。估计过程存在噪声时，估计性能下降，统计上次优化，而且计算耗时，所以一阶矩方法应用较少。

对于实际应用需要将定义式离散化，通常采用差分形式，且以中心有限差分最为常见。对于多分量信号，瞬时频率的定义就失去了它的意义。这种情况下，一般根据应用的需要，分离各个分量，再分别求出每个分量的瞬时频率。

瞬时频率提供一种作为时间函数的频域能量集中的表现形式，在信号识别、跟踪、估计和建模方面有重要意义。针对瞬时频率的时频分析方法主要从估计精度、计算复杂性、信号包含噪声的特点以及信号本身的多分量特性等方面进行研究。

2 估计方法

对于一个合理的时频分布，其能量峰脊应出现在瞬时频率的迹线附近，利用时频分布可以从两个角度估计瞬时频率：基于时间一阶条件矩的瞬时频率估计和基于谱峰检测的瞬时频率估计。基于能量检测方法进行频率估计实际上是依据式（1）进行的。已经知道这两种定义对于WVD分布是等价的，但基于一阶矩方法由于有较大的计算量，而没有额外的性能优势，所以文献中多采用能量检测的瞬时频率估计方法。

2.1 线性时频方法

常见的线性时频表示主要有短时傅里叶变换、Gabor变换以及小波变换。从考察信号的频率成分随时间的演化特性角度来说，小波变换的结果令人费解，因为它本质上是一种时间—尺度的多分辨分析。而且一旦母小波选择不当，应用效果大受影响。

短时傅里叶变换方法又称加窗傅里叶变换算法，它的基本思想是假定非平稳信号在分析窗函数的一个短的时间间隔内是准平稳，并移动分析窗函数，计算出各个不同时刻的功率谱。对频谱进行峰值检测即完成了对信号瞬时频率的估计。

该时频方法所加的时窗平滑了噪声的影响，使STFT有较好的抗噪声干扰的能力。同时该方法是线性，不会产生交叉项干扰，而且方法简单，计算量小。主要缺陷是：对于一定的时刻，只是对附近窗口内的信号作分析，存在所谓的“窗效应”，其时间和频率分辨率受窗长的约束，不能同时优化。通常将其作为粗略估计值，作为精确估计的基础。

2.2 双线性时频表示

克服时频分辨率受限的一种方法是引入双线性时频表示。这类分布多是魏格纳分布的变形，可以用统一的形式表示，习惯称之为Cohen类时频分布。

2.2.1 Wigner-Ville分布（WVD）

WVD最早由Wigner[3]于1932年提出，Ville[4]于1948年把它应用到信号处理领域。WVD定义为：

它是双线性乘积核或瞬时自相关函数关于时延的傅里叶变换。

对于常数振幅二次相位的线性调频信号，WVD方法具有理想的时频聚集性能，可以使瞬时频率曲线和时频能量峰脊重合，得到无偏的瞬时频率估计。文献[5]指出存在高斯白噪声的情况下，线性调频信号需要SNR在3 dB以上，估计性能才是优化的，并随着SNR的下降而缓慢下降；在3 dB以下，性能显著变差。

对于其他频率特征的信号，通过选择适当的窗函数，使窗内的频率特征近似线性，也可得到无偏估计的效果，但是信号必须限制在较高的信噪比水平。当信号信噪比较低时，这种方法表现出较高的方差。由此可以看出，WVD适应于线性时频信号，或者高信噪比的高阶调频信号；相应地研究WVD分布，需要克服强噪声、非线性调频以及算法固有的交叉项问题。这些问题的改善分散在下述不同的时频分析中。

WVD计算量大，没有合适的快速算法，离实时处理还有一定的距离。已有的WVD快速算法，一般采用快速傅里叶变换的方法，包括Boashash[6]等提出的重排核函数法和Martin[7]等人提出共轭对消法。

2.2.2 多项式魏格纳分布

Boashash[8]等提出了多项式魏格纳分布（PolynomialWigner-Ville Distribution，PWVD），目的在于解决多项式相位调频情况下的频率估计。PWVD定义：

其中 Kz（t，τ）为核函数：

在c=0.5，b=1，q=2的情况下，PWVD转化为WVD。

非线性调频问题转化为确定上述核函数中系数使多项式核可以将二次或多次调频信号转换成正弦波，从而使PWVD对多项式调频信号具有理想的时频聚集性。但是对于多分量信号，PWVD的交叉项将变得更加复杂，将噪声的影响在时频平面弥散化，加剧了噪声的干扰。因此，在信噪比较低时，基于PWVD的瞬时频率估计的方差较高，即使采用优化的核设计，与线性调频情况下的方差也要高出4 dB以上。另外，PWVD在实现时，必须对原信号插值，多项式核相乘使计算变得复杂。

2.2.3 交叉魏格纳分布

与WVD类似，可以定义两个信号或者一个信号的多分量的交叉魏格纳分布（Cross Wigner-Ville Distribution XWVD），即瞬时互相关函数关于时延的傅里叶变换，表达式为：

应用XWVD[9]估计瞬时频率的基本原理是在每次迭代过程这中，应用基于XWVD的谱峰检测估计出的瞬时频率来产生新的参考信号，新的参考信号和待估计的信号生成一个新的XWVD，再从新生成的XWVD估计出瞬时频率。通过迭代使信号能量聚集程度增加，从而使在噪声中估计出瞬时频率的概率增加，而且算法具有收敛性。选择合适的幅度变化率，达到C-R界的估计方差，信噪比最小可以约为-2 dB。XWVD的定义不仅克服了WVD中交叉项问题，而且提高了频率估计在低信噪比情况下的性能。缺点是不能处理多分量信号，而且是迭代算法。

2.2.4 伪 WVD（Ps.WVD）

WVD本质上是一种二次型时频分布，对多分量信号将出现交叉项，在时频平面引入了模糊，这是WVD的一个重大缺陷，主要通过核函数的设计来加以抑制。其中伪Wigner-Ville分布通过与窗函数的卷积实现平滑操作，减少交叉项影响，变换结果出现负值。虽然采用平滑伪魏格纳分布，可以得到正值结果，但是由于平滑作用使得时频面出现模糊，影响时频分辨率。另外，窗函数设计应用于STFT或者WVD，从本质上存在一个前提，即信号在窗口内满足平稳性的假设，对于快速时变信号，这个条件很难在任意时刻得到满足。因此，该类分布在时频估计方面应用较少，并出现了大量的基于自适应时频表示的频率估计方法。

2.3 自适应时频表示

2.3.1 自适应窗口长度选择

短时傅里叶变换和WVD变形中采用固定长度窗函数，很难在所有时刻满足平稳性的假设，从而引起频率估计偏差或者方差的增大，相应于大的或小的窗口长度。为了解决这个问题，可以采用自适应于估计频率的窗口长度。这种方法可以应用于SFTF或WVD的多种变形形式。基本方法一般是首先推导基于时频频谱分布进行频率估计的偏差和方差的表达式，通过均方误差的表达式在偏差和估计方差之间进行折中，从而得到自适应的窗口长度。

为解决STFT的“窗效应”问题，文献[10]中提出了一种自适应窗长的STFT方法，该窗函数为含两个控制函数的高斯窗函数，该方法对于瞬态信号和长时间信号分量有一定的自适应能力。文献[11]中提出短时线性调频窗函数的STFT对信号的瞬时频率进行估计。文献[12]利用自适应方法估计各种调频信号，通过蒙特卡罗仿真，指出除了三角形式的调频信号以外，包括阶梯状、三角函数形式的调频，其估计的累积偏差和方差都明显低于伪魏格纳分布的结果。

2.3.2 自适应核函数方法

为了更好地在抑制交叉项的同时保持很好的信号时频聚集性，对于具有不同时变频率特性的信号，核函数也要根据每个信号不同的特性来选取，每一步的最佳参数通过基函数与信号残余内积的模的最大化确定。信号自适应扩展得到基函数后，再利用不同的时频分布，利用能量峰值探测求较好的瞬时频率估计。

S.Mallat和Z.Zhang[13]提出在更大的范围内寻找匹配原始信号的基函数，利用Gauss函数的尺度变换、时频移位组成一个Gabor函数集。然后根据最大匹配投影原理寻找最佳基函数的线性组合，以达到自适应分解之目的。但是该迭代算法对时频平面的划分是一种格型分割（2j），当待分析的信号是chirp信号时，这种匹配相当于零阶逼近，势必会造成分解过程存在许多截断和分量之间的混合畸变。为克服这一缺陷，L.Angrisani等[14]提出了“小调频波变换”，简称 CT。该算法在核函数上乘以一个二次项的调频信号，由两个参数分布表示切率（chirprate）和时频平面单元的曲率，并且采用高斯函数的开方形式使其具有归一化能量。对于多项式相位的频率估计误差在1～2%以内；多分量情况下也给出了类似的结果，但是STFT和伪WVD却出现了严重的干扰，只有部分瞬时频率可以估计出来。当待分析信号的瞬时频率具有一定的周期性时，文献[15]给出与CT类似地称之为WT的时频变换。在高斯母函数的基础上，乘上一个用初始幅度、频率、相位中心频率的线性调频信号，与仿真验证该方法在处理瞬时频率具有时间周期性信号时频率估计的均方误差根的大小比CT变换要低3 dB以上。同样的信号利用STFT处理，多分量之间的相互干扰使得即使经过重排处理，其可读性也较差。

文献[16]采用高斯函数的尺度变换、旋转、时移和频移构成一簇基函数。其中旋转通过分数阶傅里叶变换来实现。通过旋转角度的调节，可使基函数在频带中心的邻域内更好地匹配信号。该方法适应于局部波形具有高斯形状的特征，具有较好的时频聚集性和抗噪声干扰的能力。但由于基于高斯函数的分数阶傅里叶变换的自适应时频分布不满足时频边缘特性，因此不能用时间一阶条件矩的方法估计瞬时频率。文献[17]在2 dB噪声的环境中，准确地估计出瞬时频率，而应用WVD估计出的瞬时频率由于受噪声的影响而出现了锯齿，特别是在两端，估计值出现振荡，且远小于实际值。

N.E.Huang提出了著名的Hilbert-Huang变换方法（简称为 HHT），其主要内容是经验模型分解 （Empirical Mode Decomposition，简称为EMD）和Hilbert谱。在这一理论中，Huang等人引入了一种单分量的数学模型，他称之为内模函数（Intrinsic Mode Function，简称为 IMF），并给出了将一个任意信号自适应地分解成IMFs之和的分解算法，即EMD。这种方法主要是利用基于经验的模式分解方法，把复杂的数据序列分解为简单的、有限个基本模式分量，得到的基本模式分量具有很好的Hilbert变换特性，使得瞬时频率具有实际的物理意义。同时对复杂数据序列的瞬时频率的引入，消除了传统的信号分析方法如傅里叶分析所产生的表示非线性，非平稳信号的伪谐波。由于它是基于信号局部特征的，因此它适用于非线性、非平稳信号的处理。

3 研究与发展趋势

从考察信号的频率成分随时间的演化特性角度来说，时频分析技术力求在计算复杂性和时频分辨率之间找到最佳平衡点。客观地说，各种时频分析技术难分轩，关键是其适合何种类型的信号。从时频分析技术应用于瞬时频率估计的研究过程中，不难发现以下的发展趋势：

1）瞬时频率与时频分布关系密切，但时频分布的计算量非常大，寻找WVD的快速计算方法，提高WVD的实现速度是一件十分有意义的跟踪，必将进一步推动时频分布在实际中的应用。

2）高阶时频分布由于较好的时频聚集性能，成为估计非线性瞬时频率的有力工具，然而对于多分量信号，高阶时频分布的交叉项将变得非常复杂。如何在尽可能保留高阶时频分布好的时频聚集性能前提下减少交叉项的干扰是一个重要的发展方向。

3）自适应时频分布的优越性能使其必将成为瞬时频率估计的趋势。如何实现自适应信号分解，以及在信号长度增加的情况下，如何实现参数优化是很有前途的研究方向。

[1]ArmstrongE.Amethodofreducingdistributesinradiosignaling by a method of frequency modulation[J].Proc.IEE， 1963（24）:669-740.

[2]Gabor D.Theory of communication[J].Proc.IEE，1946，93（3）:429-457.

[3]Qian S，CHEN Dang-pang.Joint time-frequency analysis[J].Prentice Hall，1996:41-52，181-187.

[4]Ville J.Theorie et application de al notion de signal analytique[J].Cable et Transmission，1948，2A（1）:61-74.

[5]Rao P，Taylor F J.Estimation of the instantaneous frequency using the discrete wigner distribution[J].Electronics Letters，1990， 26（4）:246-248.

[6]Boashash B.An efficient real-time implementation of the Wigner-Ville Distribution[J].IEEE Transaction on Acoustics Speech and Signal processing， 1987（35）:1611-1618.

[7]Martin W，Flandrin P.Wigner spectral analysis of nonstationary process[J].IEEE Transaction on Acoustics Speech and Signal processing， 1987（33）:1461-1470.

[8]Barkat B，Boashash B.Design of higher order polynomial wigner-ville distributions[J].IEEE transaction on signal processing， 1999，47（9）:2608-2611.

[9]Boashash B.Time-frequency signal analysis[M].In Advances in Spectrum Analysis and Array processing，1990.

[10]Czerwinski R N.Adaptive short-time Fourier analysis[J].IEEE Signal Processing Letters， 1997， 4（2）:42-45.

[11]Kwok H K，Jones D L.Improved Instantaneous frequency estimation using an adaptive short-time Fourier transform[J].IEEE Trans.Signal Processing， 2000，48（10）：2964-2972.

[12]Changdra Sekhar S，Sreenias T V.Adaptive spectrogram vs.adaptive pseudo-Wigner-Ville distribution for instantaneous frequency estimation[J].Signal Processing， 2003（83）:1529-1543.

[13]Mallat S，Zhang Z.Adaptive time-frequency transform[J].IEEE transaction on Signal processing，1993（III），241-244.

[14]Angrisani L， D’Arco M.A measurement method based on an improved version ofthe Chirplettransform for instantaneous frequency estimation[C]//IEEE Instrumentation and Measurement Technology Conference，2001:1123-1130.

[15]Angrisani L， D’Arco M， Moriello Schiano L，et al.On the use of the Warblet transform for instantaneous frequency estimation[J].IEEE，IMTC 2004：935-941.

[16]马世伟.非平稳信号的参数自适应时频表示及其应用的研究[D].上海：上海大学，2006.

[17]CHEN Guang-hua，CAO Jia-lin，WANG Jian.Use adaptive time-frequency distribution ofinstantaneous frequency estimation[J].Systems Engineering and Electronics，2002，24（1）:31-35.