孫宏軍 ，王小威

(1. 天津大學(xué)電氣自動化與信息工程學(xué)院，天津 300072；2. 天津市過程檢測與控制重點實驗室，天津 300072)

正弦信號頻率估計是機械工程、傳感器、通信及雷達(dá)等領(lǐng)域中經(jīng)常遇到的一個問題[1]．正弦信號頻率估計的關(guān)鍵是估計精度、頻率范圍和估計算法的速度．Rife等[2]提出基于最大似然(ML)算法估計正弦信號頻率，估計誤差臨近 CRLB．雖然 ML算法估計性能較好，但由于其算法復(fù)雜、速度慢，不利于實時處理，實際中一般很少直接采用 ML估計．相比而言，利用離散傅里葉變換(DFT)進(jìn)行正弦信號估計是一種更為有效的方法，得到了廣泛的應(yīng)用．DFT相當(dāng)于一組窄帶濾波器，即使在低信噪比時，只要DFT長度足夠大，仍對白噪聲有較強抑制作用．基于 DFT的譜分析方法可以利用 FFT變換實現(xiàn)快速運算，算法計算量小、速度快，特別適合實時信號處理．但由于受非同步采樣和頻率分辨率的制約[3]，基于 DFT變換的FFT算法存在頻譜泄露和柵欄效應(yīng)等不足[4]，信噪比較低時易受干擾．為了對強干擾正弦信號進(jìn)行實時、精確的頻率估計，國內(nèi)外不少學(xué)者提出了基于異頻相位處理[5]、全相位時移相位差頻率估計[6]及基于 FFT的插值算法[7]，如拋物線內(nèi)插法、雙線性對稱頻率內(nèi)插法、Rife算法等[8]．其中 Rife算法簡單、計算量小，獲得更多的關(guān)注．

1970年，Rife等[9]首先提出了利用正弦信號DFT主瓣內(nèi)兩條譜線(分別為最大值譜線和次大值譜線)的幅度變化來提高正弦信號的頻率和幅度估計精度的方法，即 Rife算法，取得了較高估計精度．Jain等[10]進(jìn)一步研究了這種算法并加以改進(jìn)，提出了插值DFT正弦信號參數(shù)估計方法，利用DFT變換主瓣內(nèi)兩條譜線的幅度比值估計信號的頻率與幅度最大值對應(yīng)頻率的相對偏差，既保留了基于 DFT的正弦信號參數(shù)估計快速的特點，又提高了正弦信號的估計精度．

Rife算法在低信噪比且實際頻率位于量化頻率點附近時，容易出現(xiàn)插值方向判斷錯誤而產(chǎn)生較大誤差．為此，Quinn[11]提出利用 DFT變換主瓣內(nèi)次大值與最大值之比的實部代替幅值之比的頻率估計方法，從而消除插值方向錯誤造成的頻率估計誤差．刁瑞朋等[12]提出適用于多頻信號的高精度頻域迭代插值(IIN)算法，先利用常規(guī)算法估算頻率參數(shù)，然后利用估算值進(jìn)行補償，再利用補償后得到的頻率參數(shù)值進(jìn)行計算，通過不斷迭代，提高估算精度和穩(wěn)定性，但該算法計算量較大．鄧振淼等[13]提出了改進(jìn)的 MRife算法，采用移頻技術(shù)進(jìn)行頻率估計，頻移量為固定值 1/3，該算法性能不隨被估計信號頻率分布而產(chǎn)生波動，在低信噪比下不存在發(fā)散問題，但有時需進(jìn)行二次移頻[14]．王宏偉等[15]提出了 I-Rife算法，利用了頻移技術(shù)和頻譜細(xì)化技術(shù)，提高了頻率修正方向的準(zhǔn)確度，但需要進(jìn)行4個點的DFT運算且仍需移頻，運算過程較復(fù)雜．孫宏軍等[16]提出了基于相角判據(jù)的 Rife算法(簡稱 P-Rife算法)，降低了方向判斷錯誤比例，整體性能優(yōu)于 Rife算法．本文在文獻(xiàn)[16]的基礎(chǔ)上提出基于幅值-相角判據(jù)的修正 Rife算法(簡稱 A-P-Rife算法)：通過設(shè)定頻移門限值，利用 FFT變換后的幅值和相角計算出頻移因子；對頻移因子滿足門限值的信號采用 Rife算法進(jìn)行估計，不滿足門限值的信號則采用相角判據(jù)進(jìn)行估計，在整個頻段上得到了穩(wěn)定且精度較高的估計性能．

1 Rife算法分析

采用未考慮高斯白噪聲背景的單一頻率復(fù)正弦信號[17]作為正弦信號簡化模型，則單頻率復(fù)正弦信號可以表示為

式中 A、fc、θc分別為正弦信號的幅度、頻率和初相角．對上述信號在[0，T]區(qū)間內(nèi)按等間隔 t＝T/N進(jìn)行采樣，T為采樣時間，得到長度為 N的序列 s(n)，則采樣得到的信號序列可表示為

序列s(n)的DFT為

S(k)幅值最大處的離散頻率索引值記作m，m＝int[fcT]，int[x]表示取最接近x的整數(shù)．DFT 的最大譜線粗測頻率為，量化頻率為，fs為采樣頻率．對于較大的 N，在幅度最大處，S(k)幅度可表示為

其中為信號頻率與其DFT幅度最大處對應(yīng)頻率的相對偏差，δ 的變化范圍為-0.5～0.5．在緊鄰m的左右兩側(cè)的兩條譜線中幅度較大處對應(yīng)的離散頻率索引值記作S(m2)的幅度可近似表示為

Am2和Am的比值記作α，當(dāng)N較大時有

根據(jù)式(6)可得|δ|的估計值[18]為

根據(jù)δ 值可對由離散頻譜得到的 fc的估計值進(jìn)行插值，從而得到更精確的頻率估計值

式(8)中符號根據(jù)m2的位置確定，若m2＝m+1則取加號，反之取減號．

可得Rife算法的正弦波頻率計算公式[19]為

在適度的信噪比條件下，當(dāng)fc位于兩相鄰頻譜的中點時與較強且幅值很接近，Rife方法頻率估計性能較好．反之，當(dāng)信噪比較低且十分接近mfs/N時的值較小，噪聲影響大，誤判概率大．當(dāng)信噪比較低趨近零時，出現(xiàn)誤判概率較大，降低了Rife算法的精度．

2 改進(jìn)Rife算法

2.1 P-Rife算法

Rife算法計算量小，當(dāng)被估計頻率位于兩相鄰量化頻譜中心區(qū)域時，能保證較高的頻率估計精度．但當(dāng)信噪比較低、趨近零時，產(chǎn)生較大頻率估計誤差[20]．P-Rife算法研究接近于零時頻譜相角的特性，從相角的角度研究正弦頻率估計．

由式(3)可知

文獻(xiàn)[16]給出的相角判據(jù)為

即檢測 k0及其左右的相角，根據(jù)式(12)確定頻移因子，并將其代入式(8)中求解頻率值．

2.2 M-Rife算法

信號頻率位于量化頻率點附近即趨近零時，Rife算法估計精度降低．修正的Rife(M-Rife)算法通過對信號進(jìn)行移頻，使新信號位于兩相鄰量化頻率點的中心區(qū)域．M-Rife算法基本思想：定義(m+1/3，m+2/3)為離散頻率點 m與 m+1之間的中心區(qū)域；先用 Rife算法進(jìn)行頻率估計得到，再判斷是否位于兩相鄰量化頻率點的中心區(qū)域；若是，則將作為最后的頻率估計值；否則需對原信號進(jìn)行移頻，使新信號的頻率位于兩個相鄰離散頻率點的中心區(qū)域，再利用Rife算法進(jìn)行頻率估計．

如果滿足式(13)時則認(rèn)為位于DFT量化頻率中心區(qū)域，并作為最終估計值；否則需要進(jìn)行修正．為使被估計信號頻率盡量接近量化頻率中點，將信號 s(n)向左或向右頻移δk量化頻率單位，δk的計算式為

當(dāng)時，r=1，此時譜線右移；反之，譜線左移．頻率計算式為

式中頻移δk是一個不確定值，與初始估計值相關(guān)．在實際操作中，取一個固定值δk＝1/3進(jìn)行移頻，計算信號頻率．

2.3 I-Rife算法

由于Am-0.5和的幅度比Am±1的幅度大，噪聲免疫力強，選取代替Am±1作為修正方向r=±1的判據(jù)更合理．利用頻移技術(shù)將信號 s(n)的頻譜進(jìn)行向左或向右移動δk量化單位，使被估計頻率盡量接近兩相鄰譜線中點，再進(jìn)行估計．頻移因子δk與修正因子Δk滿足關(guān)系Δk+δk＝1/2，則頻移因子為

I-Rife算法的頻率計算公式為

式 中 ：r=±1，當(dāng)時，r=1；當(dāng)時，r=-1．

3 A-P-Rife算法

3.1 算法思想

從前面分析可知，在Rife算法中，當(dāng)實際信號頻率 fc位于兩相鄰頻譜中點，即 fc＝fs(m+r/2)/N附近時，性能較好；當(dāng)信號頻率位于量化頻率點附近時，Rife算法估計性能較差，而相角判據(jù)估計精度較高．本文結(jié)合兩種算法性能特點，提出基于幅值-相角判據(jù)的修正Rife算法，即A-P-Rife算法．

A-P-Rife算法的基本思想為：首先設(shè)定一個頻移門限值此門限值的確定如下：根據(jù)Rife算法在被估計頻率位于兩相鄰量化頻率點的中心區(qū)域時估計性能較好，位于兩相鄰量化頻率點時且在低信噪比時估計性能較差，故所取門限值需接近兩相鄰頻譜中間區(qū)域，并結(jié)合計算機仿真實驗結(jié)果，在此門限值下，該算法估計性能穩(wěn)定且精度較高．利用頻移門限值對不同頻段采用不同的頻率估計算法，對滿足門限值的頻段采用 Rife算法進(jìn)行估計，不滿足門限值的頻譜使用相角判據(jù)，在不同頻段發(fā)揮不同算法的性能優(yōu)勢，以提高正弦頻率估計的整體性能，且計算簡便，估計精度較高．

3.2 算法實現(xiàn)

當(dāng)N足夠大時，由式(10)和式(11)可知

則 S(k)的相角由 A(k)和φ(k)確定．若δ 取值范圍為0～0.5，則在 k0處有

假定落入第一象限，則落入第三象限也落入第三象限，則有

由式(21)可知，S( k0)和S( k0-1)的相角均落入第一象限，而的相角落入第三象限，并且有

由式(22)可知S( k0)和S( k0-1)的相角相差不大，而S( k0+1)的相角有很大的跳變．當(dāng)φ(k0)落入其他象限時，也可得到類似結(jié)論．基于上述分析可以得到相角判據(jù)的頻移因子δ為

因此可知相角判據(jù)算法如下：檢測 k0及其左右的相角，并且向有跳變的方向進(jìn)行修正，即當(dāng)時修正方向r=-1；當(dāng)時修正方向r=1；當(dāng)時修正方向r=0．式(23)中絕對值符號只表示相角之間的夾角大?。嘟桥袚?jù)算法計算量較小，修正方向誤判率較低，即使在低信噪比時，也能獲得較高精度．

3.3 算法流程

A-P-Rife算法流程如圖1所示．

圖1 A-P-Rife算法流程Fig.1 Flow chart of A-P-Rife algorithm

算法步驟如下：

步驟 1對采樣序列{s(n)}進(jìn)行 FFT分析，搜索頻譜最大值位置m及左右相鄰位置m-1、m+1，并記錄各自的頻譜幅值A(chǔ)m、Am-1和Am+1；

步驟 2根據(jù)Rife算法中式(7)計算出 ?|δ|；

步驟 3設(shè)定頻移門限值并做判斷，若 ? 0.4 ≤|δ|≤0.5時，繼續(xù)步驟4，否則轉(zhuǎn)向步驟5；

步驟 4根據(jù)算出的 ?|δ|確定修正方向，并將其作為頻移因子，利用式(8)計算信號頻率；

步驟 5根據(jù) FFT分析的結(jié)果，由式(22)計算出頻譜最大值位置及左右相鄰位置處的相角值，運用相角判據(jù)算法，根據(jù)式(23)計算出頻移因子 ?|δ|及其修正方向r的值；

步驟 6利用式(8)計算信號頻率．

4 計算量分析與性能分析

4.1 算法復(fù)雜性分析

基于幅值-相角判據(jù)的修正Rife算法需要做一次N點FFT變換，還需要進(jìn)行3次反正切計算．做一次N點 FFT需要進(jìn)行(N/2)lbN次復(fù)數(shù)乘法和NlbN次復(fù)數(shù)加法．以常用的 8位單片機為例，進(jìn)行一次反正切計算需要約 12次加(減)法、5次乘法和 2次除法．考慮到當(dāng)采樣點 N較大時，計算量略大于 Rife算法，近似于P-Rife算法，但遠(yuǎn)小于M-Rife和I-Rife算法，可滿足實時性要求，如表1所示．

表1 算法計算量比較Tab.1 Comparison of computational complexity

4.2 性能分析和算法比較

許多工程現(xiàn)場環(huán)境中的熱噪聲和散粒噪聲可近似用高斯白噪聲進(jìn)行數(shù)學(xué)建模．設(shè)高斯白噪聲序列為 r(n)，服從于(0，σ,2)的高斯分布．信噪比定義為．仿真中采樣頻率 fs＝1,024,Hz，采樣點數(shù)N＝1,024，因此量化頻率f＝fs/N＝1,Hz．

對于復(fù)正弦波信號，在初相角、幅度和頻率 3個參數(shù)均未知的情況下，頻率估計標(biāo)準(zhǔn)差CRLВ定義為

基于以上分析，對正弦波加高斯白噪聲信號序列進(jìn)行 5,000次蒙特卡洛仿真實驗，以驗證基于幅值和相角判據(jù)的改進(jìn) Rife算法的性能．設(shè)f0為一量化頻率，現(xiàn)取 f0＝20,fs,/N＝20,Hz，從 f0到 f0+f,/2取 11個離散頻率點 fi＝f0+i.Δf,/20(i＝0，1，2，…，10)，對頻率fi的正弦波按 Rife、P-Rife、М-Rife、I-Rife 及本文的A-P-Rife算法進(jìn)行 5,000次 Мonte-Carlo模擬實驗.圖 2給出了 5種算法在不同信噪比下的頻率估計標(biāo)準(zhǔn)差和 CRLВ的比較．圖 3給出了 5種算法在不同信噪比下的頻率修正方向誤判率．

從圖2可以看出，其他4種算法的估計精度均優(yōu)于 Rife算法．P-Rife算法估計精度優(yōu)于 М-Rife算法，低于A-P-Rife算法和I-Rife算法，而A-P-Rife算法在估計精度上接近 I-Rife算法，逼近 CRLВ．由圖3可知，在低信噪比時，Rife算法、М-Rife算法頻率修正方向誤判率逐漸增大且大于其他3種算法.

A-P-Rife算法修正方向誤判率較小，且低于 PRife算法，修正方向準(zhǔn)確度近似于I-Rife算法．

圖2 5種算法頻率估計標(biāo)準(zhǔn)差與CRLB比較Fig.2 Comparison of five algorithms for standard deviation of frequency estimation with CRLB

圖3 5種算法頻率修正方向誤判率Fig.3 Frequency correction direction error rates of five algorithms

表2 仿真結(jié)果1(SNR＝4,dB，CRLB＝0.007,7)Tab.2 Simulation results 1(SNR＝4,dB，CRLB＝0.007,7)

表2～表4列出了5種算法的11個離散頻率點頻率估計的平均絕對誤差(МAE)、均方根誤差(RМSE)，由表中數(shù)據(jù)可看出，A-P-Rife算法在整個頻段上保持較高的估計精度與穩(wěn)定性．在較高的信噪比下，其他4種算法估計性能均優(yōu)于Rife算法，且A-P-Rife算法估計性能優(yōu)于Rife算法、P-Rife算法、М-Rife算法，接近 I-Rife算法，逼近 CRLВ；在低信噪比下，Rife算法、М-Rife算法穩(wěn)定性較低，A-PRife算法估計性能穩(wěn)定，并高于P-Rife算法，略低于I-Rife算法．

表3 仿真結(jié)果2(SNR＝0,dB，CRLB＝0.012,2)Tab.3 Simulation results 2(SNR＝0,dB，CRLB＝0.012,2)

表4 仿真結(jié)果3(SNR＝-4,dB，CRLB＝0.019,3)Tab.4 Simulation results 3(SNR＝-4,dB，CRLB＝0.019,3)

5 結(jié) 論

(1) 通過計算機模擬仿真，對5種算法進(jìn)行頻率估計，并將頻率估計標(biāo)準(zhǔn)差與CRLВ進(jìn)行比較．估計性能從高到低依次為 I-Rife、A-P-Rife、P-Rife、МRife、Rife算法，算法復(fù)雜程度由高到低依次為 IRife、М-Rife、P-Rife、A-P-Rife、Rife 算法，其中本文提出的A-P-Rife算法整體性能最佳．

(2) A-P-Rife算法降低了修正方向誤判率、減小了平均誤差、估計性能可以很好接近于CRLВ．

(3) 本文提出的 A-P-Rife算法計算量小，易于硬件實現(xiàn)，可對正弦信號實時處理，具有實際指導(dǎo)意義和應(yīng)用價值．

參考文獻(xiàn)：

［1］Kay S. A fast and accurate single frequency estimator[J]. IEEE Transactions on Acoustics，Speech，and Signal Processing，1989，37(12)：1987-1990.

［2］Rife D C，Boorstyn R R. Single tone parameter estimation from discrete-time observation[J]. IEEE Trans on Information Theroy，1974，20(5)：591-598.

［3］刁瑞朋，孟慶豐，范 虹. 基于 Rife-Vincent窗衰減信號插值離散傅里葉變換[J]. 機械工程學(xué)報，2015，51(4)：1-7.Diao Ruipeng，Meng Qingfeng，F(xiàn)an Hong. Interpolation algorithms based on Rife-Vincent window for discrete Fourier transforms of damped signals[J]. Journal of Mechanical Engineering，2015，51(4)：1-7(in Chinese).

［4］Abatzoglou T J. A fast maximum likelihood algorithm for the frequency estimation of sinusoid based on Newtons method[J]. IEEE Trans on Acoust Speech Signal Process，1985，33(1)：77-89.

［5］杜保強，周 渭. 基于異頻相位處理的高精度頻率測量系統(tǒng)[J]. 天津大學(xué)學(xué)報，2010，43(3)：262-266.Du Baoqiang，Zhou Wei. High precision frequency meas ure ment system based on different frequency phase processing[J]. Journal of Tianjin University，2010，43(3)：262-266(in Chinese).

［6］黃翔東，白瑞朋，靳旭康. 基于頻移補償?shù)娜辔粫r移相位差頻率估計[J]. 天津大學(xué)學(xué)報：自然科學(xué)與工程技術(shù)版，2017，50(6)：649-655.Huang Xiangdong，Bai Ruipeng，Jin Xukang. Frequency estimation of all phase time shifted phase difference based on frequency shift compensation[J]. Journal of Tianjin University：Science and Technology，2017，50(6)：649-655(in Chinese).

［7］許 珉，劉凌波. 基于三次樣條函數(shù)的加 Blackmanharris窗插值 FFT算法[J]. 電力自動化設(shè)備，2009，29(2)：59-63.Xu Min，Liu Lingbo. Add blackman-harris window interpolation FFT algorithm based on cubic spline function[J]. Electric Power Automation Equipment，2009，29(2)：59-63(in Chinese).

［8］王旭東，劉 渝，鄧振淼. 基于修正Rife算法的正弦波頻率估計及FPGA實現(xiàn)[J]. 系統(tǒng)工程與電子技術(shù)，2008，30(4)：621-624.Wang Xudong，Liu Yu，Deng Zhenmiao. The estimation of sine wave frequency based on the corrected Rife algorithm and FPGA implementation[J]. Journal of Systems Engineering and Electronics，2008，30(4)：621-624(in Chinese).

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

［10］Jain V K，Collins W L Jr，Davis D C. High-accuracy analog measurements via interpolated FFT[J]. IEEE Trans on Instrumentation and Measurement，1979，28(2)：113-122.

［11］Quinn B G. Estimation of frequency，amplitude and phase from the DFT of a time series[J]. IEEE Trans on Signal Processing，1997，45(3)：814-817.

［12］刁瑞朋，孟慶豐. 基于長程泄漏補償?shù)牡逯?FFT方法及應(yīng)用[J]. 振動與沖擊，2013，32(22)：1-6.Diao Ruipeng，Meng Qingfeng. Iterative interpolation FFT method and application based on long range leakage compensation[J]. Journal of Vibration and Shock，2013，32(22)：1-6(in Chinese).

［13］鄧振淼，劉 渝，王志忠. 正弦波頻率估計的修正Rife算法[J]. 數(shù)據(jù)采集與處理，2006，21(4)：473 477.Deng Zhenmiao，Liu Yu，Wang Zhizhong. Modified Rife algorithm for frequency estimation of sinusoid wave[J]. Journal of Data Acquisition and Processing，2006，21(4)：473-477(in Chinese).

［14］謝 勝，陳 航，于 平. 基于 FFT并二次修正的Rife頻率估計算法[J]. 探測與控制學(xué)報，2010，32(4)：48-53.Xie Sheng，Chen Hang，Yu Ping. Rife frequency estimation algorithm based on FFT and twice modifying[J].Journal of Detection and Control，2010，32(4)：48-53(in Chinese).

［15］王宏偉，趙國慶. 正弦波頻率估計的改進(jìn) Rife算法[J]. 信號處理，2010，26(10)：1573-1576.Wang Hongwei，Zhao Guoqing. Improved rife algorithm for frequency estimation of sinusoid wave[J]. Signal Processing，2010，26(10)：1573-1576(in Chinese).

［16］孫宏軍，徐冠群. 基于相角判據(jù)的 Rife算法的渦街信號處理方法[J]. 儀器儀表學(xué)報，2013，34(12)：2860-2866.Sun Hongjun，Xu Guanqun. Modified Rife frequency estimation algorithm based on phase criterion for vortex signal processing[J]. Chinese Journal of Scientific Instrument，2013，34(12)：2860-2866(in Chinese).

［17］齊國清，賈欣樂. 插值 FFT估計正弦信號頻率的精度分析[J]. 電子學(xué)報，2004，4(4)：625-629.Qi Guoqing，Jia Xinle. Accuracy analysis of frequency estimation of sinusoid based on interpolated FFT[J].Acta Electronica Sinica，2004，4(4)：625-629(in Chinese).

［18］胥嘉佳，劉 渝，鄧振淼，等. 正弦波信號頻率估計快速高精度遞推算法的研究[J]. 電子與信息學(xué)報，2009，31(4)：865-869.Xu Jiajia，Liu Yu，Deng Zhenmiao，et al. The study of frequency estimation algorithm of the sin signal with fast speed and high accuracy[J]. Journal of Electronics and Information Technology，2009，31(4)：865-869(in Chinese).

［19］周喜慶，趙國慶，王 偉. 實時準(zhǔn)確正弦波頻率估計綜合算法[J]. 西安電子科技大學(xué)學(xué)報：自然科學(xué)版，2004，31(5)：657-660.Zhou Xiqing，Zhao Guoqing，Wang Wei. Real-time and accurate frequency estimation algorithm of the sine wave signal[J]. Journal of Xidian University：Natural Science，2004，31(5)：657-660(in Chinese).

［20］Xiao Y C，Wei P，Xiao X C，et al. Fast and accurate single frequency estimator[J]. Electronics Letters，2004，40(14)：910-911.