復合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