复合FFT插值*

陈奎孚， 赵建柱

（1.中国农业大学理学院74# 北京，100083） （2.中国农业大学工学院 北京，100083）



复合FFT插值*

陈奎孚1， 赵建柱2

（1.中国农业大学理学院74# 北京，100083） （2.中国农业大学工学院 北京，100083）

现有的插值FFT一般仅利用谱峰附近最高/次高谱线对。为了降低估计方差，首先利用谱峰附近4条谱线给出3个估计式，然后对其加权平均。给出了使方差最小的优化权系数，以及优化后的方差。理论分析表明：新方法的方差小于Quinn方法；在整周期采样的情形下，新方法方差只有Quinn方法的4/9。采用仿真算例对新方法进行了考核，结果表明：在高信噪比（SNR，50 dB）且非整周期采样条件下，模拟方差与理论解一致；就所模拟的范围（SNR低至0 dB），模拟方差超过理论值最大仅有25％。

参数估计；频谱；快速傅里叶变换（FFT）；方差；优化

引 言

尽管已经发展了很多形形色色的参数估计方法，但基于傅里叶变换的谱方法仍是处理周期信号的重要的方法。其中一个重要原因就是存在快速傅里叶变换（fast Fourier transform，简称FFT）。然而直接由FFT谱线读出的信号参数的误差有时很大，在极端情形下，幅值偏差可达31％，初相位偏差可达90°。这种显著的偏差在FFT刚提出不久就被发现了。随后的发展就是长期努力不懈地构造各式各样的窗函数来抑制上述偏差［1-2］。

另外一条思路是利用谱峰附近谱线采用类似参数模型的方法算出参数。Rife［3］针对矩形窗给出了利用幅值比的计算公式，并被推广到其他窗函数［4］。笔者发现只有对Rife-Vincent窗函数才存在简单的计算公式［5］，但为了控制估计方差，在通信领域和电子领域一般不加窗。然而不加窗又会出现插值方向错误［6-7］。此外，Rife计算公式在整周期采样条件下方差比半周期采样的方差大。为了降低估计方差，刘渝教授［8］及其同事给出基于迭代的修正Rife算法，并对迭代的初值进行了细致研究［9］，最近给出了递推算法［10］。Quinn［11］则用第3条谱线的相位的方法来降低方差，笔者也给出了3条谱线法［7］。由于笔者给出的方法基于复比值而非幅值［12-13］（当然还有基于相位的方法［14］和基于实部的方法［15］），所以很容易推广到多谱线［5］。

本研究探索如何利用4条谱线来提高FFT插值的统计性能。

1 复单频信号频谱模型

笔者仅考虑如下的复单频模型（实单频信号信号可以看作为两个复单频信号之和）：

其中：xk（k=0～N-1）为采样序列；N为采样点数，也就是FFT的长度；Δt为采样间隔；A，α和φ分别为所感兴趣成分的幅值、角频率和初相位；zk（k=0～N-1）为零均值复白噪声序列。

物理上不存在复噪声，但是用希尔伯特变换将实信号转成解析信号时，就会在数学上出现复噪声。形式上其中实部和虚部满足如下的白噪声条件：

加窗就是对记录的时间段内信号沿时间轴乘以不同权重，以增强化谱分析的期望特性；而矩形窗相当于不加窗，也就是直接使用原始数据做谱分析，笔者仅考

为了进行分析，笔者假定噪声在概率意义上相对于信号幅度A很小。这样可将式（13）用泰勒公式对噪声展开，保留到一阶小量为

从式（14）可以看出，无论l为何值（靠近主瓣），式（12）都为α的无偏估计。

2 复合平均的方差

通常插值只用最高/次高谱线对（图1中k+1～k+2）。非整周期采样在实际操作中难以避免，这样频谱的泄露使得除了最高/次高谱线外，还会在其周围出现其它较明显的谱线。考察图1中的k～k+3四条谱线，相应地可以有3个估计式，即

图1　谱峰附近示意图Fig.1 Profile around the spectrum peak

3 优化方差

图2显示了βL，min和βR，min随δ变化的趋势。对该图有如下评述：

图2　权系数随频率偏移的变化Fig.2 Dependence of optimal weighting coefficients on the frequency shift

当δ=0，中间估计式αM权系数最大，左右两个αL和αR的贡献很小。应指出，此时βL，min和βR，min均为很小的负值（-1/82），相应地αM权系数为42/ 41。

随着δ自中心向两边偏移，两侧估计式的权系数加大，如向右偏移，βR，min正向增大，αM权系数减小。当δ=1/2时，βL，min=0。这是因为谱线正好落在k+2条谱线上，而用于估计αL的两条谱线k和k+1完全为噪声（因为整周期采样使得这两条谱线上没有任何能量），因而其贡献自然应为零。还应指出的是，当δ=1/2时βR，min=4/9≠1/2。但这时离散谱以k+2谱线对称，如果只考虑αR和αM两者平均的话，直觉上应该二者各占一半。然而βR和βL是被同步优化的，因此βL趋近于零的极限和βL先验地等于零不同。正因为βL趋近于零的干扰对αR和αM程度不同，才使得αR和αM的优化权系数有差异。

当δ从1/2变化到1时，αR权系数持续增加，而αM权系数进一步下降。当δ=1时，信号频率正落在k+2和k+3两条谱线中间，αR精度最高，相应权系数也最大。但一般来说，这种情形出现可能性比较小。因为通常的最高/次高谱线的选择默认了真实频率在k+1和k+2谱线之间，即|δ|＜0.5。

当δ＞0时，αL权系数很小；反之当δ＜0，αR权系数也很小。

对某些δ，权系数可能小于零，而统计学经常使用的权系数为正值，比如有文献用幅值加权的方法提高精度［17］，但能否简单地应用于提高频率精度需要进一步研究。因为这里理论上的优化权函数可能是正数，也可能是负数。

优化的权系数随δ而变化，而后者是需要估计的量，因此必须找到计算权系数的近似方法。其中一种方法是利用最高/次高谱线估计出一个近似值，然后根据式（29）和（30）计算权系数。在图2中还画出了两条曲线

它们与βL，min在感兴趣的区间上|δ|＜0.5比较接近，因此可以用来近似βL，min。在实际操作中需要用取样值Yk～Yk+3代替式中相应的Xk～Xk+3。βR，min取βL，min的镜象对称即可。

由式（28）可以得到所关心的最小方差

图3给出了新方法与3条谱线法［7］，以及Quinn方法［11，19］的比较。从图可以看出，在整周期采样条件下（δ=1/2）下，新方法的误差显著低于Quinn方法。后者的而新方法为前者的4/9。

如同Quinn估计，当δ=1/2时新方法的性能最差，这对应整周期采样条件；反之新方法在半周期采样条件（δ=0）下的性能最好。当稍低于Quinn的这是因为后者已经非常接近方差的下限Cramer-Rao界了。

图3　最小方差的理论解Fig.3 The theoretical minimum variance for compared estimators

4 仿真比较和分析

本节将比较3种方法：Quinn方法［11，19］，简单复比值法［12］和新给出的复合加权平均法。按照式（1）生成仿真信号，采样间隔Δt=1 s，序列长度N= 256，这样Δω=2π/（NΔt）=0.024 5 rad/s。频率α从34.5Δω等间隔扫描到35.5Δω，间隔为0.025Δω。由于初相位对估计方差没有影响，所以每个样本序列的初相位在［0，2π］均匀分布（通过对MATLAB的rand函数平移和扩张来实现）。共计考察了6个噪声级别，即SNR=50，25，10，5，3和0 dB。每个信噪比×频率共计模拟20 000个序列。

生成仿真信号后直接进行FFT①笔者未对生成的信号去零。显然仿真序列的均值未必精确等于零，如果对它强制去零，那么虚假的均值会导致谱泄露。在信噪比很高情况下，估计的方差将由虚假均值造成的谱泄露所主导。，按照图1所示选择谱线，分别用3种方法估计。新方法权系数按式（26）计算，其中的δ根据最高/次高谱线对确定。对20 000个序列分别统计3种方法的模拟方差最后的结果如图4所示。

3种方法中，简单复比值法的方差最大。在SNR比较高的图4（a）中除δ=1/2外，简单复比值法与Quinn理论表达式一致。但随SNR降低，不吻合的区域自δ=1/2向两侧对称扩展。然而，它却远远低于简单幅值比法中因插值方向错误所造成的方差［19］。

Quinn方法的模拟方差与其理论表达式基本吻合。图4（a）的δ=0和δ=1处模拟方差比理论值稍大。对SNR比较低且δ＞0.5的情形，模拟方差比理论值略大，其原因需进一步研究。

图4　模拟方差与理论值的比较Fig.4 Comparisons of the empirical variance against theoretical predictions

新方法表现最好。SNR最高的图4（a）表明除了δ=1/2外，其模拟方差与理论表达式吻合良好。随着SNR降低，模拟方差大于理论值的区域自中心δ=1/2向外扩展，但前者最大也只有后者的125％（就所模拟的数据）。差异的原因可能是由于分析所采用的泰勒展开的截断误差，该误差相对于噪声的影响随SNR减小而上升。应该指出的是：用于复合的3个估计均相当于简单复比值法，而中间估计式方差就对应图4中的空心小圆；3者经过优化平均之后，方差显著降低。

5 结束语

为了提高信号估计精度，笔者针对矩形窗，研究了利用谱峰附近4条谱线来估计参数的复合FFT插值方法。连续4条谱线给出了3个估计式，新方法就是对其的加权平均。在理论上导出了使方差最小的权系数，并用仿真序列对理论进行了检验。

理论分析表明新方法的方差小于Quinn方法，降低误差效率与偏离周期采样的程度有关。效率最高发生于整周期采样，此时新方法方差只有Quinn方法的4/9。优化的权函数随采样偏离整周期程度而变化。笔者提供了两组经验权系数。

采用仿真方法比较了Quinn方法、简单复比值法和新给出方法。考核结果显示，在高信噪比（SNR =50 d B）且非整周期采样条件下，新方法的仿真方差与理论值基本吻合，而就所模拟的范围（SNR低至0 d B），模拟方差超过理论值最大仅有25％。

［1］ Harris F J.On the use of windows for harmonic analysis with the discrete Fourier transform［J］.Proc IEEE，1978，66（1）：51-83.

［2］ Rejin I S，Reljin B D，Papic V D.Extremely flat-top windows for harmonic analysis［J］.IEEE Transaction on Instrument and Measurement，2007，56（3）：1025-1041.

［3］ 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.

［4］ Offelli C，Petri D.Weighting effect on the discrete time Fourier transform of noisy signals［J］.IEEETransactions on Instrumentation and Measurement，1991，40（6）：972-981.

［5］ Chen K F，Jiang J T，Crowsen S.Against the longrange spectral leakage of the cosine window family［J］. Computer Physics Communications，2009，180（6）：904-911.

［6］ 齐国清，贾欣乐.插值FFT估计正弦信号频率的精度分析［J］.电子学报，2004，32（4）：625-629. Qi Guoqing，Jia Xinle.Accuracy analysis of frequency estimation of sinusoid based on interpolated FFT［J］. Acta Electronica Sinica，2004，32（4）：625-629.（in Chinese）

［7］ Yang X Z，Li H Y，Chen K F.Optimally weighted average of the interpolated fast fourier transform in both directions［J］.IET Science，Measurement and Technology，2009，3（2）：137-147.

［8］ 邓振淼，刘渝，王志忠.正弦波频率估计的修正Rife算法［J］.数据采集与处理，2006，21（4）：474-477. Deng Zhenmiao，Liu Yu，Wang Zhizhong.Modified rife algorithm for frequency estimation of sinusoid wave［J］.Journal of Data Acquisition&Processing，2006，21（4）：474-477.（in Chinese）

［9］ 邓振淼，刘渝.正弦波频率估计的牛顿迭代方法初始值研究［J］.电子学报，2007，35（1）：104-107. Deng Zhenmiao，Liu Yu.The starting point problem of sinusoid frequency estimation based on newton′s method［J］.Acta Electronica Sinica，2007，35（1）：104-107.（in Chinese）

［10］胥嘉佳，刘渝，邓振淼.正弦波信号频率估计快速高精度递推算法的研究［J］.电子与信息学报，2009，31（4）：865-869. Xu Jiajia，Liu Yu，Deng Zhenmiao.A research of fast and accurate recursive algorithm for frequency estimation of sinusoid signal［J］.Journal of Electronics&Information Technology，2009，31（4）：865-869.（in Chinese）

［11］Quinn B G.Estimating frequency by interpolation using Fourier coefficients［J］.IEEE Transactions on Signal Processing，1994，42（5）：1264-1268.

［12］陈奎孚，王建立，张森文.频谱校正的复比值法［J］.振动工程学报，2008，21（3）：314-318. Chen Kuifu，Wang Jianli，Zhang Senwen.Spectrum correction based on the complex ratio of discrete spectrum around the main-lobe［J］.Journal of Vibration Engineering，2008，21（3）：314-318.（in Chinese）

［13］Li Y F，Chen K F.Eliminating the picket fence effect of the fast fourier transform［J］.Computer Physics Communications，2008，178（7）：486-491.

［14］黄翔东，王兆华.基于全相位频谱分析的相位差频谱校正法［J］.电子与信息学报，2008，30（2）：293-297. Huang Xiangdong，Wang Zhaohua.Phase difference correcting spectrum method based on all-phase spectrum analysis［J］.Journal of Electronics&Information Technology，2008，30（2）：293-297.（in Chinese）

［15］黄玉春，黄载禄，黄本雄.基于FFT滑动平均极大似然法的正弦信号频率估计［J］.电子与信息学报，2008，30（4）：831-835. Huang Yuchun，Huang Zailu，Huang Benxiong.FFTBased moving average maximum likelihood single-tone frequency estimation［J］.Journal of Electronics&Information Technology，2008，30（4）：831-835.（in Chinese）

［16］Schoukens J，Renneboog J.Modeling the noise influence on the Fourier coefficients after a discrete Fourier transform［J］.IEEE Transactions on Instrumentation and Measurement，1986，IM-35（3）：278-286.

［17］庞浩，李东霞，俎云霄.应用FFT进行电力系统谐波分析的改进算法［J］.中国电机工程学报，2003，23（6）：50-54. Pang Hao，Li Dongxia，Zu Yunxiao.An improved algorithm for harmonic analysis of power system using FFT technique［J］.Proceedings of the CSEE，2003，23（6）：50-54.（in Chinese）

［18］Rife D C，Boolstyn R R.Single-tone parameter estimation from discrete-time observation［J］.IEEE Transactions on Information Theory，1974，20（5）：591-598.

［19］ 齐国清.几种基于FFT的频率估计方法精度分析［J］.振动工程学报，2006，19（1）：92-96. QI Guoqing.Accuracy analysis and comparison of some FFT-based frequency estimators［J］.Journal of Vibration Engineering，2006，19（1）：92-96.（in Chinese）

TN911.6

10.16450/j.cnki.issn.1004-6801.2015.02.000

陈奎孚，男，1969年12月生，博士、教授。主要研究方向为振动工程和生物物理。曾发表（《Against the long-range spectral leakage of the cosine window family，computer physics communication》2009年第180卷第6期）等论文。

E-mail：Chen KuiFu@Hotmail.com

简介：赵建柱，男，1963年1月生，副教授。主要研究方向为从事车辆动力学。

E-mail：zhjzh@cau.edu.cn.

2014-01-01；

2014-04-28