打开文本图片集
摘要：为了进一步提高加性高斯白噪声背景中正弦信号的频率估计精度，提出了一种新的基于插值快速傅里叶变换（FFT）的正弦信号频率估计算法。首先，对N点正弦采样序列进行等长度时域补零延长，再进行2N点FFT； 然后， 搜索幅度最大离散谱线位置得到频率粗估计值； 最后， 采用幅度最大谱线以及原信号的离散时间傅里叶变换（DTFT）在幅度最大谱线左右两侧的两点抽样值进行精估计。仿真结果表明，当信号实际频率位于FFT两条离散谱线之间任意位置时，所提算法的频率估计均方根误差均接近克拉美罗下限，具有较好的一致性，估计精度高于Candan算法、Fang算法、三谱线合理结合（RCTSL）算法和Aboutanios算法， 且信噪比阈值较低，估计性能优于现有频率估计算法。
关键词：频率估计；快速傅里叶变换；正弦信号；离散时间傅里叶变换；数字信号处理
中图分类号：TN911.72
文献标志码：A
0引言
加性高斯白噪声背景中的正弦信号频率估计是数字信号处理领域中的一个经典课题，目前在通信、雷达、声呐以及电子对抗等领域得到了广泛应用。频率估计一般考虑3个因素：频率估计精度、估计范围以及算法实现所需的运算量。文献[1]中提出的最大似然估计法估计误差达到克拉美罗下限（CramerRao Lower Bound，CRLB），是最优估计估计算法，但是算法复杂，运算量大，无法实时处理。基于快速傅里叶变换（Fast Fourier Transform，FFT）的频率估计算法具有运算速度快、对正弦信号有显著的信噪比增益等优点，适合对速度要求高的实时处理系统，所以此类估计算法引起了国内外学者的广泛关注[2-11]。
由于FFT得到的是离散频率值，可以将信号频率表示为f0=（km+δ）Δf。其中： km是FFT幅度最大谱线对应的离散频率索引值； δ表示信号频率与FFT幅度最大谱线位置的相对频率偏差，δ∈[-0.5，0.5]； Δf是FFT的频率分辨率。基于插值FFT的正弦信号频率估计主要分为两个步骤，粗估计和精估计。其中粗估计通过搜索FFT最大谱线完成，精估计借助一定的插值算法得到信号真实频率与粗估计之间的相对频率偏差，各种插值FFT算法的差异仅体现在精估计中。在正弦信号幅度和频率都未知的情况下，利用其FFT离散频谱主瓣内幅度最大值和次大值可以解出正弦信号的频率（和幅度）[2-5]。直接利用FFT主瓣内两条谱线的幅度估计信号频率算法简单、运算量小。但是，当信号实际频率与FFT幅度最大离散谱线对应的频率接近（δ接近0）时，由于FFT主瓣内次大谱线幅度很小、容易受噪声影响导致频率估计误差增加[6]。因此，又出现了多种对其进行改进的方法[7-12]。文献[7]中Jacobsen等提出利用FFT频谱最大的三根谱线对频率估计值进行校正，在低信噪比和FFT点数较少时结果较好，但估计精度不高。在文献[8]中Candan对Jacobsen算法的系数进行了修正，提高了估计精度，但通过仿真实验发现，其估计误差依然较大，不适合对精度要求较高的场合。在文献[9]中Candan对文献[8]的算法进行了改进消除了估计偏差，改善了原估计算法在高信噪比时的性能，但其在低信噪比时的频率估计均方误差仍明显高于CRLB。Fang等在文献[10]中提出，对N点正弦采样序列进行N点时域补零后，采用2N点FFT频谱中最大谱线相邻的两根谱线幅度对相对频率偏差进行估计。在时域对原信号进行补零延长后，FFT的离散谱线之间的间隔变小，该方法利用2N点FFT离散频谱中最大值左右两侧的谱线幅度估计信号频率，由于信号时域补零延长后频谱主瓣变宽，FFT在主瓣内的谱线数增加，频谱最大值左右两侧（离散频率值+/-1）的谱线幅度比原来（信号未延长）的N点FFT主瓣内的第二大谱线更靠近最大谱线、幅度更大，不容易受噪声影响，因而频率估计精度有所提高。但该算法在相对频率偏差的整个取值范围内估计性能不稳定。文献[11]提出了三谱线合理结合（Rational Combination of Three Spectrum Lines，RCTSL）算法，对N点正弦采样序列进行N点时域补零后，根据最小二乘法估计思想，采用2N点FFT频谱中最大的三根谱线幅度对相对频率偏差进行估计，估计精度高，算法复杂度较低。但RCTSL算法在推导过程中，用到的能量谱函数的一阶泰勒级数展开，省去了高阶项，导致估计误差的产生。离散信号x（n）的FFT的离散频谱X（k）只出现在离散频率k取整数值的位置。文献[12]将其推广到离散频率k取非整数值的情况，直接利用FFT离散频谱最大值左右离散频率值+0.5处和-0.5处的频谱幅度来求解信号频率。其基本思想也是尽量避免利用幅度很小的FFT系数，提高抗噪声干扰的性能，并通过多次迭代以降低频率估计误差。该方法是现有估计算法中精度最高的，由于每次迭代需要多次计算FFT系数，增加了计算量。受文献[10]和[12]的启发，本文提出一种新的正弦信号频率估计算法，在对原信号补零加长1倍进行2N点FFT找到幅度最大离散谱线位置后，直接计算原信号的离散时间傅里叶变换（DiscreteTime Fourier Transform，DTFT）在幅度最大离散谱线位置+0.5Δf和-0.5Δf处的抽样值，即X0.5与X-0.5，并根据X0、X0.5和X-0.5得到信号频率的估计值。由于这里的X0.5和X-0.5比文献[12]中的X0.5和X-0.5更靠近原信号的DTFT的幅度最大值，因此受噪声影响更小，在噪声背景中频率估计性能比前述方法更好。
第11期
樊磊等：基于快速傅里叶变换的正弦信号频率高精度估计算法
计算机应用 第35卷
1本文算法
在不考虑噪声时，单一频率复正弦信号可表示为：
x（t）=Aej（2πft+θ0）（1）
式中：A、 f、θ0分别为复正弦信号的幅度、频率和初相。
经过采样率为fs，采样点数为N的采样后，信号的离散化形式可以表示为：
x（n）=A exp[j（2π（f/fs）n+θ0）]；n=0，1，2，…，N-1（2
对x（n）进行N点时域补零延长得到x′（n），并进行2N点FFT，得到：
X[k]=∑N-1n=0Aejφej2πffsne-jπNnk=
Aejφejπ（N-1）（ffs-k2N）sin （πN（f/ fs - k/2N））sin （π（f/ fs-k/2N））（3
X[k]幅度最大值处的离散频率索引值为km，根据X（k）幅度最大值处的位置km可以得到信号频率的粗估计值为kmΔf， Δf=fs/（2N）为FFT的频率分辨率，即相邻谱线间的间隔。
为了描述方便，将2N点FFT系数表示为：
Xp=X[km+p]=∑2N-1n=0x′（n）e-jπnkm+pN（4）
其中：p表示与谱线峰值位置km的间隔。时域补零相当于频域插值，所以FFT频谱的谱线间隔减小为原来的1/2，此时的谱线间隔代表的实际频率为Δf=fs/（2N）。而FFT频谱的主瓣宽度为2fs/N，为谱线间隔的4倍，FFT频谱主瓣内出现4根谱线X0、X1、X-1和X2（或者是X-2，视δ的情况而定）。
将f0=（km+δ）fs/（2N）代入式（3），并进行整理可以得到：
Xp=aejθ02N1-ejπ（δ-p）1-ejπ（δ-p）/N（5
虽然幅度最大离散谱线X0及相邻两根谱线X1、X-1幅度较大，且均包含了相对频率偏差δ的信息，可以用来对δ进行估计。但是在δ=0.5附近时，X1与X-1之一幅度比较小，较容易受到噪声影响，导致估计误差较大。虽然离散信号x′（n）的2N点FFT的离散频谱X[k]只出现在离散频率k取整数值的位置，但为了进一步减小估计误差，可以考虑采用幅度最大离散谱线X0以及x′（n）的DTFT在X0位置+0.5Δf和-0.5Δf处的抽样值来对δ进行估计。x′（n）的DTFT在X0位置+0.5Δf和-0.5Δf处的抽样值表示如下
X（ejω）"f=（km±0.5）Δf=DTFT[x′（n）]|f=（km±0.5）Δf=
∑2N-1n=0x′（n）e-j2πfn|f=（km±0.5）Δf=X±0.5（6
根据DTFT与FFT的关系，上式中的两个抽样值，就是按照式（4）计算的X0.5和X-0.5。如图1所示，在无噪声情况下，x′（n）的DTFT幅度谱为图中的实线包络，其峰值对应频率为信号频率f，幅度最大离散谱线X0对应离散频率索引值为km（图中km=4）。从图中可以看出，X0.5和X-0.5比X1和X-1距离幅度最大离散谱线X0更近，幅度也更大，更不容易受到噪声影响。
图片
图1无噪声时x′（n）的DTFT幅度谱与2N点FFT系数
（N=8，A=1，δ=0.2， f0=（N/4+δ）fs/N
根据文献[12]当N的取值远大于π（δ-p）时，ejπ（δ-p）/N≈1+jπ（δ-p）/N，式（5）可表示为：
Xp=bpδ-p（7）
其中
bp=aejθ02π（1-ejπ（δ-p））（8
分别计算b-0.5、b0.5和b0如下：
b-0.5=aejθ02π（1-ejπ（δ+0.5））=aejθ02π（1-jejπδ）
b0.5=aejθ02π（1-ejπ（δ-0.5））=aejθ02π（1+jejπδ）
b0=aejθ02π（1-ejπδ）
于是有：
X-0.5=b-0.5δ+0.5=aejθ02π（1-jejπδ）δ+0.5（9
X0.5=b0.5δ-0.5=aejθ02π（1+jejπδ）δ-0.5（10
X0=b0δ=aejθ02π（1-ejπδ）δ（11
整理后可以得到：
X0.5X0=δδ-0.5·1+jejπδ1-ejπδ（12
X-0.5X0=δδ+0.5·1-jejπδ1-ejπδ（13
对式（12）、式（13）进行整理得到：
X0.5X01-0.5δ-1=ejπδ[X0.5X01-0.5δ+j]（14
X-0.5X01+0.5δ-1=ejπδ[X-0.5X01+0.5δ-j]（15
上面两式相除，并进行整理得到
X0.5（δ-0.5）-X0δX-0.5（δ+0.5）-X0δ=X0.5（δ-0.5）+jX0δX-0.5（δ+0.5）-jX0δ（16
可以得到δ的估计表达式
δ^=0.5[（1-j）X0.5+（1+j）X-0.5]（1-j）X0.5-（1+j）X-0.5+2jX0（17
由于噪声的影响，式（17）右边的表达式的计算结果可能存在虚部。为了保证得到的δ^是实数，对式（17）右边的表达式取实部，得到：
δ^=Real0.5[（1-j）X0.5+（1+j）X-0.5]（1-j）X0.5-（1+j）X-0.5+2jX0（18
于是得到信号频率的估计表达式：
f^=（km+δ^）Δf（19
为了进一步提高估计精度，可以采用本文算法对上述频率估计结果进行修正。首先采用本文算法对信号频率进行估计，得到相对频率偏差δ的初步估计值δ^1，令km′=km+δ^1。此时km′Δf接近信号实际频率f，Xkm′+0.5与Xkm′-0.5的幅度均较大，且二者幅度很接近，不容易受到噪声影响。接着在式（4）中用km′取代km，采用本文算法得到δ的估计值δ^2，用δ^2对初步估计值δ^1进行修正。则信号频率的估计表达式如下：
f^=（km+δ^1+δ^2）Δf（20
上述算法步骤如下：
1）对N点的采样序列x（n）进行N点时域补零延长；
2）进行2N点FFT，得到X[k]；
3）搜索最大谱线位置，将kmΔf作为信号频率的粗估计值；
4）采用本文算法，根据式（18）得到相对频率偏差δ的初步估计值δ^1；
5）按式（4）计算Xδ^1、Xδ^1+0.5与Xδ^1-0.5；
6）采用本文算法，将Xδ^1、Xδ^1+0.5与Xδ^1-0.5代入式（18），得到δ^1的修正估计值 δ^2；
7）根据式（20）得到信号频率最终估计值。
2仿真实验及性能分析
为了验证上面提出的新的频率估计算法的性能，本章通过计算机Monte Carlo模拟实验进行仿真，并与Candan算法[8]、Fang算法[10]、 RCTSL算法[11] 以及Aboutanios算法[12]进行比较。在不补零的FFT插值算法中，Candan算法的性能是最优的[8]，所以只选择与这一种不补零的FFT插值算法进行比较。这几种算法对δ的估计表达式如下所示。
Candan算法[8]：
δ^=
tan （π/N）π/NReal{X（km-1）-X（km+1）2X（km）-X（km-1）-X（km+1）}
Fang算法[10]：
δ^=tan （π/N）π/N|X（km+1）|-|X（km-1）||X（km+1）|+|X（km-1）|
RCTSL算法[11]：
δ^=
Nπ|X（km+1）|2-|X（km-1）|2u（|X（km+1）|2+|X（km-1）|2）+v|X（km）|2；
u=64N（π5+32π），v=uπ4
Aboutanios迭代算法[12]：
设定δ估计初值δ^0=0， i=1∶Q，迭代Q次
Xp=∑N-1n=0x（n）exp{-j2πN（k0+δ^i-1+p）n}；
p=±0.5
δ^i=δ^i+12Real{X0.5+X-0.5X0.5-X-0.5}
加性高斯白噪声背景中的单一频率复正弦信号为：
x（t）=Aej（2πft+θ0）+w（t） （21）
其中w（t）是方差为σ2的高斯复白噪声，其均值为0。信噪比定义为SNR=A2/σ2。对于复正弦信号，在幅度、频率、相位3个参数都未知的情况下，频率估计的克拉美罗方差下限[1]为：
min var（f^0）=32π2T2N（N2-1）·SNR （22
当N=512，SNR=3dB， f0=（N/4+δ）Δf时，本文算法及前述几种算法频率估计的均方根误差（Root Mean Square Error，RMSE）与CRLB的比值随相对频率偏差δ的变化情况如图2所示，其中δ为进行时域补零延长后的相对频率偏差。从图中可以看出，Fang算法[10]与RCTSL算法[11]在δ的整个取值范围内估计性能不稳定，当在δ=0.5附近时，由于2N点FFT频谱中最大谱线相邻的两根谱线之一幅度比较小，较容易受到噪声影响，导致估计误差较大（对于RCTSL算法，δ=0.5时除外）。除了当δ=0.5时，本文算法频率估计的RMSE在δ的整个取值范围内比其他算法更近接近CRLB，且性能稳定。Aboutanios迭代算法[12]的频率估计精度在其他算法中是最高的，本文算法的估计精度高于Aboutanios迭代算法。
图片
图2几种算法的RMSE随δ的变化情况
当N=128， f0=（N/4+δ）fs/N，δ均匀分布在[-0.5，0.5]上时，本文算法及前述几种算法的RMSE随SNR的变化情况如图3和图4所示。从图中可以看出，Candan算法[8]的频率估计精度明显低于其他算法。本文算法的信噪比阈值与Fang算法以及RCTSL算法相同，约为-7dB，低于Candan算法和Aboutanios迭代算法。
图片
图3几种算法的RMSE随信噪比变化情况
图片
图4几种算法的RMSE随信噪比变化情况
3结语
本文提出了一种新的基于插值FFT的正弦信号频率估计算法，利用补零加长FFT频谱中的幅度最大谱线以及原信号的DTFT在幅度最大谱线左右两侧的两点抽样值得到正弦信号频率的估计公式。仿真结果表明，该算法的频率估计性能优于几种最新的基于FFT幅度插值的频率估计算法。
与Candan算法相比，本文算法的频率估计RMSE明显更低，信噪比阈值也更低。与Fang算法以及RCTSL算法相比，本文算法的性能具有较好的一致性。不论实际信号频率位于FFT两条离散谱线之间的什么位置，频率估计RMSE都接近CRLB，而Fang算法和RCTSL算法当δ接近0.5时，其RMSE明显高于CRLB（对于RCTSL算法，δ=0.5时除外）。与Aboutanios迭代算法相比，本文算法采用的X0.5和X-0.5比Aboutanios迭代算法中的X0.5和X-0.5更靠近原信号的DTFT的幅度最大值，受噪声影响更小，因此频率估计RMSE更低，信噪比阈值也更低。虽然补零延长后进行FFT以及计算原信号的DTFT的抽样值增加了一些运算量，但带来的好处是更高的估计精度和更好的抗噪性能。该算法适用于对估计精度要求较高的场合。
参考文献：
[1]RIFE D C， BOORSTYN R R. Singletone parameter estimation from discretetime observations[J]. IEEE Transactions on Information Theory， 1974， 120（5）： 591-598.
[2]RIFE D C， VINCENT G A. Use of the discrete Fourier transform in the measurement of frequencies and levels of tones[J]. Bell System Technical Journal， 1970， 49（2）： 197-228.
[3]JAIN V K， COLLINS W L， DAVIS D C. Highaccuracy analog measurements via interpolated FFT[J]. IEEE Transactions on Instrumentation and Measurement， 1979， 28（2）： 113-122.
[4]QUINN B G. Estimating frequency by interpolation using Fourier coefficients[J]. IEEE Transactions on Signal Processing， 1994， 42（5）： 1264-1268.
[5]QUINN B G. Estimation of frequency， amplitude and phase from the DFT of a time series [J]. IEEE Transactions on Signal Processing， 1997， 45（3）： 814-817.
[6]QI G， JIA X. Accuracy analysis of frequency estimation of sinusoid based on interpolated FFT[J]. Acta Electronica Sinica， 2004，32（4）： 625-629. （齐国清，贾欣乐. 插值FFT估计正弦信号频率的精度分析[J]. 电子学报，2004，32（4）： 625-629.）
[7]JACOBSEN E， KOOTSOOKOS P. Fast， accurate frequency estimators [J]. IEEE Signal Processing Magazine， 2007， 24（3）： 123-125.
[8]CANDAN C. A method for fine resolution frequency estimation from three DFT samples [J]. IEEE Signal Processing Letters， 2011， 18（6）： 351-354.
[9]CANDAN C. Analysis and further improvement of fine resolution frequency estimation method from three DFT samples [J]. IEEE Signal Processing Letters， 2013， 20（9）： 913-916.
[10]FANG L， DUAN D， YANG L. A new DFTbased frequency estimator for singletone complex sinusoidal signals [C]// MILCOM 2012： Proceedings of the 2012 IEEE Military Communications Conference. Piscataway： IEEE， 2012：1-6.
[11]YANG C， WEI G. A noniterative frequency estimator with rational combination of three spectrum lines[J]. IEEE Transactions on Signal Processing， 2011， 59（10）： 5065-5070.
[12]ABOUTANIOS E， MULGREW B. Iterative frequency estimation by interpolation on Fourier coefficients [J]. IEEE Transactions on Signal Processing， 2005， 53（4）： 1237-1242.
相关热词搜索： 正弦 算法 变换 频率 信号