复白噪声中复正弦波频率估计方法研究

王 乐 苏小敏 杜 林 李春化

(西安电子工程研究所 西安 710100)

1 引言

对于随机信号，由于其无始无终、能量无限，对其作傅里叶变换不收敛，因而不能像确定性信号那样确定此类信号的频谱。虽然此类信号能量无限，但是其功率未必无限，因此对于随机信号而言，常用功率谱描述其频率特性。功率谱估计分为两大类:参数化方法和非参数化方法。非参数化方法依赖于传统的傅里叶变换法，其包括 BT法、周期图法等［1］，参数化方法包括MUSIC算法、AR模型算法、最大熵等估计算法。利用MUSIC算法、AR模型、最大熵等现代谱估计算法可以对正弦信号频率进行精确的估计［2］，但是这些算法复杂、计算量大，难以实时处理从而影响其应用。

对于任意形式的确定性信号x(n)，对其进行傅里叶分解，都可以分解为直流与许多余弦(或正弦)分量线性组合［1］。然而在基于离散傅里叶变换(DFT)的经典谱估计中，由时域截断引起的频谱泄露是不可避免的，能量泄露会引起较大的估计误差［3］，这时我们需要对频谱进行校正以期获得精确的谱峰位置。本文提出一种功率加权求平均的频谱校正算法，并通过仿真实验在基于FFT分析基础上与单谱线法、Rife法进行对比，然后讨论了观察时间对FFT分析的影响以及信号采样频率的选择，最后针对短持续时间信号，给出频谱细化方法(Zoom-FFT)，并通过实验与FFT方法进行对比。

2 频谱校正算法

设正弦信号为:

式中，A、f0、φ0代表信号的振幅、频率和初相。对其采样，得到离散序列:

其中，Δt=1/fs代表采样间隔，fs为采样频率，M为采样点数。对(2)式进行离散傅里叶变换(DFT)得:

直接计算DFT，其计算量是与变换区间长度N的平方成正比的。当N较大时，计算量非常的大，因此直接用DFT算法进行谱分析和信号的实时处理是不切实际的。工程实际中常用快速傅里叶变换(FFT)来计算DFT，其离散频谱谱峰值对应的频率作为信号频率的初步估计，即:

其中，k表示频谱峰值对应的位置。由于其仅仅利用了峰值谱线，又被称为单谱线幅度法。

我们知道，对一个信号作DFT，前提要求其必须是有限长度的，因此必须对信号作加窗处理，窗函数的引入会导致频谱泄漏。又因为作离散傅里叶变换相当于对整个频率范围(0～fs)内作N点均匀采样，在均匀采样过程中我们并不一定能采得幅度峰值，这就导致了k取值的不准确性。该频率估计方法分辨率较低，最大频率误差可以达到 ± fs/2N［2］。

对于基于基带FFT法利用单根谱线估计频率造成误差过大问题，人们提出了诸多频谱校正方法。其中较著名、效果较好的是双线幅度法(Rife法)。Rife法利用最大谱线和次大谱线对单谱线方法进行修正，当信噪比较高时估计性能很好。但当信噪比较低或者待估计频率靠近最大谱线时，会导致次大谱线判断的失误，造成估计误差比仅仅用单谱线粗略估计误差还要大，此时则要利用插值细化技术对Rife法作进一步改进，关于 Rife法详见参考文献［4～5］。

本文提出了一种利用最大谱线及其邻近谱线，对其进行功率加权并取平均的方法来弥补单谱线估计的不足。方法示意图见图1。

图1 功率加权求平均校正频谱

图1 中虚线部分表示频谱包络，其中频谱峰值对应频率为f0，假设在其左右各取两根相邻谱线，分别对应频率值为 f-2、f-1、f1、f2，所取五根谱线如图所示，对应幅值分别为 l-2、l-1、l0、l1、l2，则信号频率估值为:

上述主要介绍了在整个频域(0～fs)范围内，对序列信号进行频域上一维搜索，找出频谱峰值，用其所对应的频率值作为原信号的频率估计值。这种估计方法原理清晰、易于理解，但计算量较大。若要提高估计精度，延长观测时间是一种行之有效的方法。采用MATLAB仿真工具进行验证，假设采样频率fs=2048Hz，采样间隔 Ts=1/fs，信号频率 f0=125.65Hz，信号幅度A=2.0，噪声设为正态分布随机噪声，方差为0.4，信噪比SNR=10dB，当观测时间分别为 t1=0.1s、t2=0.5s、t3=1s、t4=1.5s、t5=3s，采用单谱线法、Rife法和功率加权求平均方法所估计频率为，结果表 1 所示。

表1 频谱校正对比实验估计频率结果

估计频率值减去频率真值，所得差值再取绝对值，得到各方法的估计误差，如图2所示。

图2 校正频谱后的估计误差

由图2可以看出:a.双线幅度法(Rife法)并不一定优于单谱线幅度法;

b.当增加观察时间时，栅栏效应引起的频谱泄漏影响减小，频率估计误差减小，估计频率精度增高。

当观察时间固定，采样频率提高时，采样点数也相应增多，但这种情况下，频率分辨率未增加，频率估计效果是不会提高的。这是因为，虽然采样频率提高会带来采样点数N的增加，但同样导致频率轴上的折叠频率(fs/2)增加，频率分辨率fs/N未变。实验验证，假设采样频率 fs1=400Hz、fs2=800Hz、fs3=1200Hz、fs4=1600Hz、fs5=2000Hz，采样间隔 Ti=1/fsi，信号频率 f0=125.65Hz，信号幅度 A=2.0，噪声设为正态分布随机噪声，方差为0.4，信噪比SNR=10dB，当采样时间t=1s，采用上述单谱线法、Rife法和功率加权求平均方法所估计频率为、，结果如表2所示。

表2 频率递增对比实验估计频率结果

估计频率值减去频率真值，所得差值再取绝对值得到各方法的估计误差，如图3所示。

由图3可以看出，采样频率变化时频率估计误差变化不超过0.5Hz，这种些微变化可以认为是由随机噪声引起的，而与采样频率递增无关。

然而通过上面仿真实验，可以看出采样频率提高，无论哪种估计方法误差变化趋势是一样的，这就提示我们频率估计误差跟采样频率之间并不是毫无关联的。

在频率估计中，采样率的选择实际上是有一定限制的，不能随意更改采样频率的大小，否则会造成估计精度变差。换句话说，就是固定的采样频率只适合对一定频率范围内的信号进行采样，否则会给频率估计带来问题，以下进行仿真验证。

图3 频率递增后的估计误差

假设采样频率fs=100MHz，采样间隔Ts=1/fs，信号幅度A=1.0、噪声Ψ(n)为正态分布随机噪声，方差为0.01，信噪比SNR=20dB，信号观测时间t=1μs，信号频率为f，则信号可以表示为:

通过仿真，ω取100个不同值所对应的频率估计实验结果如图4所示，图中横坐标表示ω在0～π之间取值，纵坐标表示频率估计误差均方根，图4(右)为去掉图4(左)中始末两点所得。

图4 频率估计误差均方根随信号频率变化

由图4可以看出:当信号角频率ω位于0.2π～0.8π之时，即信号频率时，频率估计误差均方根较小，估计效果较好，当信号频率时，频率估计均方根非常大，估计不稳定，这验证了采样频率的选取是有一定限制的结论，因此我们在作频率估计时，涉及到采样频率时应该使2.5f＜fs＜10f，其中fs的上下限并不严格限定为10f和2.5f，仅代表fs在取值上注意脱离危险区域。

以上论述了几种较简单的频率校正方法，通过仿真实验证明了延长观测时间有利于频率估计结果的结论，并且讨论采样频率与信号频率之间的关系。

然而用增加观察时间来提高频率估计精度的方法并不是放之四海而皆准的，首先观察时间的增加势必带来计算量的增大;其次在短持续时间信号中，观测时间是有限的，不能通过增加观测时间来提高频率估计精度，这种情况下就需要改进算法，才能获得理想的频率估计精度。

3 频谱细化算法

在许多实际应用中，我们并不是对整个频谱感兴趣，往往只关注其中的一个窄频带，并且需要对这一窄频带作局部放大以进行细致的观察，为此我们应该在所关心的频带内增加谱线，以增大谱线密度［6］。近年来，频谱细化技术得到迅猛发展，常见的方法有基于复调制的Zoom-FFT法、线性调频Z变换法(Chirp-Z变换，简称CZT)、Yip-Zoom法、相位补偿细化法等。但从分析精度、计算效率、分辨率、灵活性等方面来看，Zoom-FFT法是一种非常有效的方法［7］，在工程中也得到了非常广泛的应用。

Zoom-FFT法(以下简称ZFFT法)是提供高分辨率FFT分析的一种技术，它允许在给定的一组FFT谱线中选择起始和终止频率，在所得的窄频带上能以指定的、足够高的采样频率去分析，可以成倍的提高频率分辨率，这样就可以对我们感兴趣的频率段进行细化。

在实际应用中，我们通常对信号做基带FFT估计，得到信号的粗测频率，然后以此粗测频率作为窄频率带中心频率，向左右延展确定起始频率和终止频率，起始频率和终止频率之间的窄频带作为细化分析频带，确定细化倍数，移频再进行复FFT处理。Zoom-FFT处理过程流程见图5所示。

图5 Zoom-FFT处理框图

算法具体步骤可以归纳为:

a.复调制移频，即将需要细化的频带中心(FFT粗测所得)移至频率轴原点，通常我们用数字下变频来实现;

b.低通滤波，我们需要滤出感兴趣的频段信号，加以分析。假设频率细化倍数为d，则低通滤波器的截止频率fc=fs/2d;

c.抽样。信号被调制移频和低通滤波之后，所得到的信号频带变窄，因而可以用较低的采样频率f's=fs/d进行重采样，即对原采样点每隔固定点抽样一次;

d.复FFT处理和频率回调。抽样之后得到复序列信号，对其进行复FFT处理，此时分辨率提高了d倍，然后进行频率调整，将谱线移至实际频率处即可得到细化的频谱。

仿真实验假设采样频率fs=100MHz，采样间隔Ts=1/fs，信号幅度A=2.0、噪声为正态分布随机噪声，方差为0.4，信噪比SNR=10dB，信号采样时间t=1μs，信号频率 f=30.123MHz，频带中心 fa=30MHz，细化倍数为 d=10，Zoom-FFT分析带宽为25MHz～35MHz，用基带FFT和ZFFT法进行频率估计，仿真结果见图5和图6所示。图中横坐标表示频率轴，纵坐标表示幅值。用功率加权求平均的方法进行频谱校正，做100次Monte-carlo实验得到仿真结果:FFT方法估计结果=30.215MHz，同真值比较，误差为Δ=0.092MHz;ZFFT方法频率估计为，误差为

由FFT与ZFFT频率估计对比实验可以看出，ZFFT方法频率估计在信号主瓣内采样点数明显多于FFT方法，这些意味着频谱分辨率更高，频率估计精度更大。

4 总结

频谱校正算法可以很好的修正由获取最大谱线位置不准确所带来的频率估计误差。本文的主要工作是:首先，通过实验对基于FFT分析的单谱线、双谱线和功率加权等频谱校正方法进行对比，得出功率加权算法优于单谱线法和双谱线法，这是因为功率加权法利用了更多的谱线信息;其次，通过实验验证了延长观测时间利于频率估计以及采样频率与信号频率之间的关系;最后，针对短持续时间信号观测时间有限的特点，介绍了一种行之有效的频谱细化方法，并通过仿真实验与FFT方法作对比，证明了ZFFT算法估计效果优于FFT算法。

［1］陆光华，彭学愚，张林让，毛用才.随机信号处理［M］.西安:西安电子科技大学出版社，2003.

［2］王洪先，王玉军，王宏伟.宽带自相关接收机中的一种频率估计方法［J］.火控雷达技术，2011，40(1):52-54.

［3］李天昀，葛临东.基于ZFFT和Chirp-Z的频谱细化分析中能量泄露的研究［J］.电子对抗技术，2003，18(5):11-15.

［4］STEVE KAY.A Fast and Accurate Single Frequency Estimator［J］.IEEE TRANS.ON ACOUSTICS，1989，37(12):1987-1989.

［5］王纪强，欧攀，张春熹，林志立.基于频偏校正的频率估计算法误差分析［J］.北京航空航天大学学报，2010，36(7):849-852.

［6］王力，张冰，徐伟.基于MATLAB复调制ZOOMFFT算法的分析和实现［J］.舰船电子工程，2006，4:119-121.

［7］丁康，潘成灏，李巍华.ZFFT与Chirp-Z变换细化选带的频谱分析对比［J］.振动与冲击，2006，25(6):9-12.