復(fù)合FFT插值*

陳奎孚， 趙建柱

（1.中國農(nóng)業(yè)大學(xué)理學(xué)院74# 北京，100083） （2.中國農(nóng)業(yè)大學(xué)工學(xué)院 北京，100083）

?

復(fù)合FFT插值*

陳奎孚1， 趙建柱2

（1.中國農(nóng)業(yè)大學(xué)理學(xué)院74# 北京，100083） （2.中國農(nóng)業(yè)大學(xué)工學(xué)院 北京，100083）

現(xiàn)有的插值FFT一般僅利用譜峰附近最高/次高譜線對。為了降低估計方差，首先利用譜峰附近4條譜線給出3個估計式，然后對其加權(quán)平均。給出了使方差最小的優(yōu)化權(quán)系數(shù)，以及優(yōu)化后的方差。理論分析表明：新方法的方差小于Quinn方法；在整周期采樣的情形下，新方法方差只有Quinn方法的4/9。采用仿真算例對新方法進(jìn)行了考核，結(jié)果表明：在高信噪比（SNR，50 dB）且非整周期采樣條件下，模擬方差與理論解一致；就所模擬的范圍（SNR低至0 dB），模擬方差超過理論值最大僅有25％。

參數(shù)估計；頻譜；快速傅里葉變換（FFT）；方差；優(yōu)化

引 言

盡管已經(jīng)發(fā)展了很多形形色色的參數(shù)估計方法，但基于傅里葉變換的譜方法仍是處理周期信號的重要的方法。其中一個重要原因就是存在快速傅里葉變換（fast Fourier transform，簡稱FFT）。然而直接由FFT譜線讀出的信號參數(shù)的誤差有時很大，在極端情形下，幅值偏差可達(dá)31％，初相位偏差可達(dá)90°。這種顯著的偏差在FFT剛提出不久就被發(fā)現(xiàn)了。隨后的發(fā)展就是長期努力不懈地構(gòu)造各式各樣的窗函數(shù)來抑制上述偏差［1-2］。

另外一條思路是利用譜峰附近譜線采用類似參數(shù)模型的方法算出參數(shù)。Rife［3］針對矩形窗給出了利用幅值比的計算公式，并被推廣到其他窗函數(shù)［4］。筆者發(fā)現(xiàn)只有對Rife-Vincent窗函數(shù)才存在簡單的計算公式［5］，但為了控制估計方差，在通信領(lǐng)域和電子領(lǐng)域一般不加窗。然而不加窗又會出現(xiàn)插值方向錯誤［6-7］。此外，Rife計算公式在整周期采樣條件下方差比半周期采樣的方差大。為了降低估計方差，劉渝教授［8］及其同事給出基于迭代的修正Rife算法，并對迭代的初值進(jìn)行了細(xì)致研究［9］，最近給出了遞推算法［10］。Quinn［11］則用第3條譜線的相位的方法來降低方差，筆者也給出了3條譜線法［7］。由于筆者給出的方法基于復(fù)比值而非幅值［12-13］（當(dāng)然還有基于相位的方法［14］和基于實(shí)部的方法［15］），所以很容易推廣到多譜線［5］。

本研究探索如何利用4條譜線來提高FFT插值的統(tǒng)計性能。

1 復(fù)單頻信號頻譜模型

筆者僅考慮如下的復(fù)單頻模型（實(shí)單頻信號信號可以看作為兩個復(fù)單頻信號之和）：

其中：xk（k=0～N-1）為采樣序列；N為采樣點(diǎn)數(shù)，也就是FFT的長度；Δt為采樣間隔；A，α和φ分別為所感興趣成分的幅值、角頻率和初相位；zk（k=0～N-1）為零均值復(fù)白噪聲序列。

物理上不存在復(fù)噪聲，但是用希爾伯特變換將實(shí)信號轉(zhuǎn)成解析信號時，就會在數(shù)學(xué)上出現(xiàn)復(fù)噪聲。形式上其中實(shí)部和虛部滿足如下的白噪聲條件：

加窗就是對記錄的時間段內(nèi)信號沿時間軸乘以不同權(quán)重，以增強(qiáng)化譜分析的期望特性；而矩形窗相當(dāng)于不加窗，也就是直接使用原始數(shù)據(jù)做譜分析，筆者僅考

為了進(jìn)行分析，筆者假定噪聲在概率意義上相對于信號幅度A很小。這樣可將式（13）用泰勒公式對噪聲展開，保留到一階小量為

從式（14）可以看出，無論l為何值（靠近主瓣），式（12）都為α的無偏估計。

2 復(fù)合平均的方差

通常插值只用最高/次高譜線對（圖1中k+1～k+2）。非整周期采樣在實(shí)際操作中難以避免，這樣頻譜的泄露使得除了最高/次高譜線外，還會在其周圍出現(xiàn)其它較明顯的譜線?？疾靾D1中的k～k+3四條譜線，相應(yīng)地可以有3個估計式，即

圖1　譜峰附近示意圖Fig.1 Profile around the spectrum peak

3 優(yōu)化方差

圖2顯示了βL，min和βR，min隨δ變化的趨勢。對該圖有如下評述：

圖2　權(quán)系數(shù)隨頻率偏移的變化Fig.2 Dependence of optimal weighting coefficients on the frequency shift

當(dāng)δ=0，中間估計式αM權(quán)系數(shù)最大，左右兩個αL和αR的貢獻(xiàn)很小。應(yīng)指出，此時βL，min和βR，min均為很小的負(fù)值（-1/82），相應(yīng)地αM權(quán)系數(shù)為42/ 41。

隨著δ自中心向兩邊偏移，兩側(cè)估計式的權(quán)系數(shù)加大，如向右偏移，βR，min正向增大，αM權(quán)系數(shù)減小。當(dāng)δ=1/2時，βL，min=0。這是因?yàn)樽V線正好落在k+2條譜線上，而用于估計αL的兩條譜線k和k+1完全為噪聲（因?yàn)檎芷诓蓸邮沟眠@兩條譜線上沒有任何能量），因而其貢獻(xiàn)自然應(yīng)為零。還應(yīng)指出的是，當(dāng)δ=1/2時βR，min=4/9≠1/2。但這時離散譜以k+2譜線對稱，如果只考慮αR和αM兩者平均的話，直覺上應(yīng)該二者各占一半。然而βR和βL是被同步優(yōu)化的，因此βL趨近于零的極限和βL先驗(yàn)地等于零不同。正因?yàn)棣翷趨近于零的干擾對αR和αM程度不同，才使得αR和αM的優(yōu)化權(quán)系數(shù)有差異。

當(dāng)δ從1/2變化到1時，αR權(quán)系數(shù)持續(xù)增加，而αM權(quán)系數(shù)進(jìn)一步下降。當(dāng)δ=1時，信號頻率正落在k+2和k+3兩條譜線中間，αR精度最高，相應(yīng)權(quán)系數(shù)也最大。但一般來說，這種情形出現(xiàn)可能性比較小。因?yàn)橥ǔ５淖罡?次高譜線的選擇默認(rèn)了真實(shí)頻率在k+1和k+2譜線之間，即|δ|＜0.5。

當(dāng)δ＞0時，αL權(quán)系數(shù)很??；反之當(dāng)δ＜0，αR權(quán)系數(shù)也很小。

對某些δ，權(quán)系數(shù)可能小于零，而統(tǒng)計學(xué)經(jīng)常使用的權(quán)系數(shù)為正值，比如有文獻(xiàn)用幅值加權(quán)的方法提高精度［17］，但能否簡單地應(yīng)用于提高頻率精度需要進(jìn)一步研究。因?yàn)檫@里理論上的優(yōu)化權(quán)函數(shù)可能是正數(shù)，也可能是負(fù)數(shù)。

優(yōu)化的權(quán)系數(shù)隨δ而變化，而后者是需要估計的量，因此必須找到計算權(quán)系數(shù)的近似方法。其中一種方法是利用最高/次高譜線估計出一個近似值，然后根據(jù)式（29）和（30）計算權(quán)系數(shù)。在圖2中還畫出了兩條曲線

它們與βL，min在感興趣的區(qū)間上|δ|＜0.5比較接近，因此可以用來近似βL，min。在實(shí)際操作中需要用取樣值Yk～Yk+3代替式中相應(yīng)的Xk～Xk+3。βR，min取βL，min的鏡象對稱即可。

由式（28）可以得到所關(guān)心的最小方差

圖3給出了新方法與3條譜線法［7］，以及Quinn方法［11，19］的比較。從圖可以看出，在整周期采樣條件下（δ=1/2）下，新方法的誤差顯著低于Quinn方法。后者的而新方法為前者的4/9。

如同Quinn估計，當(dāng)δ=1/2時新方法的性能最差，這對應(yīng)整周期采樣條件；反之新方法在半周期采樣條件（δ=0）下的性能最好。當(dāng)稍低于Quinn的這是因?yàn)楹笳咭呀?jīng)非常接近方差的下限Cramer-Rao界了。

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

4 仿真比較和分析

本節(jié)將比較3種方法：Quinn方法［11，19］，簡單復(fù)比值法［12］和新給出的復(fù)合加權(quán)平均法。按照式（1）生成仿真信號，采樣間隔Δt=1 s，序列長度N= 256，這樣Δω=2π/（NΔt）=0.024 5 rad/s。頻率α從34.5Δω等間隔掃描到35.5Δω，間隔為0.025Δω。由于初相位對估計方差沒有影響，所以每個樣本序列的初相位在［0，2π］均勻分布（通過對MATLAB的rand函數(shù)平移和擴(kuò)張來實(shí)現(xiàn)）。共計考察了6個噪聲級別，即SNR=50，25，10，5，3和0 dB。每個信噪比×頻率共計模擬20 000個序列。

生成仿真信號后直接進(jìn)行FFT①筆者未對生成的信號去零。顯然仿真序列的均值未必精確等于零，如果對它強(qiáng)制去零，那么虛假的均值會導(dǎo)致譜泄露。在信噪比很高情況下，估計的方差將由虛假均值造成的譜泄露所主導(dǎo)。，按照圖1所示選擇譜線，分別用3種方法估計。新方法權(quán)系數(shù)按式（26）計算，其中的δ根據(jù)最高/次高譜線對確定。對20 000個序列分別統(tǒng)計3種方法的模擬方差最后的結(jié)果如圖4所示。

3種方法中，簡單復(fù)比值法的方差最大。在SNR比較高的圖4（a）中除δ=1/2外，簡單復(fù)比值法與Quinn理論表達(dá)式一致。但隨SNR降低，不吻合的區(qū)域自δ=1/2向兩側(cè)對稱擴(kuò)展。然而，它卻遠(yuǎn)遠(yuǎn)低于簡單幅值比法中因插值方向錯誤所造成的方差［19］。

Quinn方法的模擬方差與其理論表達(dá)式基本吻合。圖4（a）的δ=0和δ=1處模擬方差比理論值稍大。對SNR比較低且δ＞0.5的情形，模擬方差比理論值略大，其原因需進(jìn)一步研究。

圖4　模擬方差與理論值的比較Fig.4 Comparisons of the empirical variance against theoretical predictions

新方法表現(xiàn)最好。SNR最高的圖4（a）表明除了δ=1/2外，其模擬方差與理論表達(dá)式吻合良好。隨著SNR降低，模擬方差大于理論值的區(qū)域自中心δ=1/2向外擴(kuò)展，但前者最大也只有后者的125％（就所模擬的數(shù)據(jù)）。差異的原因可能是由于分析所采用的泰勒展開的截斷誤差，該誤差相對于噪聲的影響隨SNR減小而上升。應(yīng)該指出的是：用于復(fù)合的3個估計均相當(dāng)于簡單復(fù)比值法，而中間估計式方差就對應(yīng)圖4中的空心小圓；3者經(jīng)過優(yōu)化平均之后，方差顯著降低。

5 結(jié)束語

為了提高信號估計精度，筆者針對矩形窗，研究了利用譜峰附近4條譜線來估計參數(shù)的復(fù)合FFT插值方法。連續(xù)4條譜線給出了3個估計式，新方法就是對其的加權(quán)平均。在理論上導(dǎo)出了使方差最小的權(quán)系數(shù)，并用仿真序列對理論進(jìn)行了檢驗(yàn)。

理論分析表明新方法的方差小于Quinn方法，降低誤差效率與偏離周期采樣的程度有關(guān)。效率最高發(fā)生于整周期采樣，此時新方法方差只有Quinn方法的4/9。優(yōu)化的權(quán)函數(shù)隨采樣偏離整周期程度而變化。筆者提供了兩組經(jīng)驗(yàn)權(quán)系數(shù)。

采用仿真方法比較了Quinn方法、簡單復(fù)比值法和新給出方法?？己私Y(jié)果顯示，在高信噪比（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］.電子學(xué)報，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］.數(shù)據(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］.電子學(xué)報，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］.電子與信息學(xué)報，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］陳奎孚，王建立，張森文.頻譜校正的復(fù)比值法［J］.振動工程學(xué)報，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］.電子與信息學(xué)報，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］.電子與信息學(xué)報，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］龐浩，李東霞，俎云霄.應(yīng)用FFT進(jìn)行電力系統(tǒng)諧波分析的改進(jìn)算法［J］.中國電機(jī)工程學(xué)報，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］.振動工程學(xué)報，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月生，博士、教授。主要研究方向?yàn)檎駝庸こ毯蜕镂锢?。曾發(fā)表（《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月生，副教授。主要研究方向?yàn)閺氖萝囕v動力學(xué)。

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

2014-01-01；

2014-04-28