基于改進的Rife測頻算法及其FPGA實現*

司偉建，郝 鑫，趙忠凱，陳 濤

（哈爾濱工程大學信息與通信工程學院，哈爾濱 150001）

0 引言

正弦波頻率估計是信號參數估計中的經典問題，也是電子對抗領域研究的重點課題之一。文獻[1]提出了利用DFT進行頻率估計的同時在時域加類sinc窗函數的頻譜校正方法，這種方法加寬了窗譜寬度，使得兩條相鄰DFT譜線產生測頻模糊，因此需要DFT點數較大才可以保證測頻精度。文獻[2]利用對相位長短延時的方法，短的時間延遲可以用于消除頻率模糊，長的延遲用于獲得較高的頻率分辨率的方法來提高測頻精度，但是這種方法受噪聲影響較大，在頻帶邊緣處測頻誤差較大。文獻[3]給出了最大似然頻率估計算法，該算法能夠達到頻率估計的下限：克拉美－羅限（CRB），但該算法計算量大，不利于工程實現。文中在分析Rife算法的基礎上，從FPGA實現的角度出發，提出了一種改進的Rife算法，仿真和實測結果均表明該算法能在短數據條件下實現對信號頻率的快速、精確測量，具有較高的工程應用價值。

1 Rife算法分析

文獻[4－5]采用最大譜線及其相鄰的次大譜線進行插值，用插值的結果進行頻率估計，稱這種算法為Rife算法。Rife算法的原理如下：

設接收的信號并被采樣后為：s（n）＝x（n）＋w（n），其中x（n）為接收到的純信號，w（n）為混入信號中的高斯白噪聲序列，其均值為0，方差為σ2。設x（n）＝aej（2πfc／fs＋φ0），其中a為信號幅度；fc為接收信號頻率；φ0為信號的初相；fs為信號采樣頻率。將采樣后的信號s（n）做N點FFT運算，最大譜線值記為｜X（k0）｜，次大譜線值記為｜X（k0＋r）｜。則Rife算法所得到的頻率估計式為：

圖1 Rife算法示意圖

其中，當｜X（k0＋1）｜≥｜X（k0－1）｜時，r＝1；相反，當｜X（k0＋1）｜＜｜X（k0－1）｜時，r＝－1。圖1為r＝1時的情況。 由于0≤所以插值的偏移量在[k0－1／2，k0＋1／2]范圍內。將偏移量與r的值代入式（1）即可估計出輸入信號的頻率值。

2 改進Rife算法分析

Rife算法通過判斷FFT峰值兩側譜線幅度的大小來確定插值的方向，通過計算來選取偏移量大小。由此可見利用FPGA實現時，需要進行除法運算，而除法運算需要占用較多的FPGA邏輯單元，運算周期長，限制了測頻速度。

經過仿真發現：將輸入信號加窗后，FFT 最大值兩側的兩個次大值的功率差G＝｜G（k0＋1）｜－｜G（k0－1）｜（其中G代表功率差）與在[k0－1／2，k0＋1／2]偏移量保持線性關系（其中k0為FFT幅度最大值所對應的下標索引值），如圖2所示。因此只要求出任意一段的線性斜率a，并且知道FFT的最大值下標索引值k0，利用功率差值G就可以由式（2）確定出輸入信號的頻率。

圖2 FFT下標索引值與功率差值關系圖

為了實現上述算法，首先需要求出斜率a的大小，在具體實現中可以利用Matlab計算出其結果，利用其結果便可以在硬件中實現改進Rife算法測頻。具體步驟如下：

1）對采樣序列進行加窗，加窗的長度與做FFT的點數相同；

2）加窗后的序列做FFT運算，取幅度峰值｜X（k0）｜兩側的兩個值｜X（k0＋1）｜和｜X（k0－1）｜；

3）將｜X（k0＋1）｜和｜X（k0－1）｜轉換成對應的功率值，并求出功率差：G＝｜G（k0＋1）｜－｜G（k0－1）｜，記錄G與k0；

4）在待測輸入信號頻率區間內，輸入信號每次按照相同的步進頻率重復1）～3），直至結束；

5）選取任意一段處于[k0－1／2，k0＋1／2]之間的頻率部分，根據3）記錄的G與k0的關系做出圖2，并求出該段斜率的倒數，即a值的大小；

6）利用5）求出的a和FFT的幅度峰值的下標索引值k0，以及功率差G，代入式（2）便可以求出輸入信號的頻率值。

3 仿真分析

設置仿真參數為：輸入信號頻率范圍：fc＝227.5＋Δfc，Δfc取[－fs／2 N，fs／2 N]頻率間隔為0.1MHz，采樣頻率為240MHz，采樣點數為64點，信號初始相位φ0取[0，2π]隨機數，噪聲均值為0，方差σ2＝0.5。對每個頻率點做1000次Monte Caro仿真實驗，當輸入信號為實信號時，對應的克拉美－羅限（CRB）由式（3）確定[7]。

圖3為SNR＝6dB時Rife算法和改進Rife算法的測頻誤差仿真結果，可以看出在｜Δfc｜≤1／3時，Rife算法的測頻誤差較大，最大均方誤差為0.7MHz，而此時改進的Rife算法的測頻誤差較小。相反當｜Δfc｜＞1／3時，改進Rife算法的測頻誤差變大，最大均方誤差為0.5MHz。

圖4為不同輸入信號信噪比條件下，改進Rife算法和Rife算法的測頻性能仿真結果。其中SNR變化范圍為：0～20dB，每個信噪比下做1000次Monte－Caro試驗，每次試驗信號輸入頻率為fc＝277.5＋Δfc（Δfc取[－fs／2 N，fs／2 N]內的隨機數）。由圖可以看出在信噪比變大時，Rife算法與改進Rife算法的測頻精度都在提高，但改進Rife算法的測頻誤差始終小于Rife算法。

圖3 SNR＝6dB時Rife算法和改進Rife算法測頻誤差比較

圖4 不同輸入信噪比下Rife算法與改進Rife算法測頻誤差比較

4 改進Rife算法的FPGA實現

將改進的Rife算法在Altera公司的EP3SL200F1152I4FPGA中進行了實現。選用的窗函數為海明窗，長度為64點，量化后存入ROM中。利用計數器產生讀ROM地址，并在每個時鐘上升沿將ROM中所對應的窗函數順次讀出。在Quartus II中調用乘法器來完成輸入信號與窗函數的乘積運算。結果送給FFT IP核，FFT點數為64點，輸出的結果轉換成幅度并查表轉換成功率，取出功率峰值處的下標索引k0，并將索引值為k0－1和k0＋1處的功率值調用減法器做差，得到功率的差值。由功率差值做地址在FPGA中進行查表求出在一個區間內的頻率偏移量，同時利用功率峰值下標索引k0及斜率a就可以計算出輸入信號的頻率。

除了主要模塊部分外，還包括控制FFT模塊進數時序以及延時控制等部分。本次選取窗函數量化位數為11位，FFT模塊輸入方式為streaming模式。

4.1 實測數據與仿真對比結果

信號源采用Agilent E4438C，利用噪聲源NC6110A對輸入信號加入相應功率的噪聲，ADC采樣速度為240MHz。輸入信號頻率為277.5 MHz，輸入信噪比SNR分別為0dB、3dB、6dB、9dB，脈沖寬度為2μs，對脈沖前沿樣64點進行頻率估計。利用signaltapII采樣4K輸出結果，利用Matlab軟件計算其均方誤差，實驗結果如表1所示。由表1可以看出，實測均方誤差稍大于理論均方誤差，其原因是：對頻率進行量化會引入量化誤差，同時ADC采樣不理想也會引入一定的誤差，對實驗結果產生一定的影響。

表1 實測結果

4.2 FPGA資源消耗

在Quartus II中建立工程，選取FPGA芯片為EP3SL200F1152I4，經過綜合并布線后消耗總資源為4%，各部分資源消耗情況如表2所示。

表2 改進的Rife算法資源消耗表

經過實測，單獨編譯一個64點FFT IP核所消耗的資源為2%，可見改進的Rife算法占用資源很少。

5 結論

針對Rife算法測量誤差大、硬件實現速度慢的特點，提出了一種改進的Rife測頻算法。從工程應用的角度，對算法FPGA實現的關鍵步驟進行了詳細闡述。該算法能夠實現對低信噪比、短數據條件下的未知信號進行快速、準確測頻，在0dB時測頻誤差小于0.6MHz／σ，處理速度可達240MHz，適合應用于電子對抗、偵察等領域。

[1]王曉君，陳禾，仲順安.改進的脈沖雷達信號加窗DFT頻譜校正算法[J].北京理工大學學報，2008 ，28（2）：164－167.

[2]趙亮，程少云.寬帶數字接收機測頻方法的研究[J].微計算機信息，2008，24（19）：292－294.

[3]Rife D C，Boorstyn R R.Single－tone parameter estimation from discrete time observations[J].IEEE Trans.Inform.Theory，1974，20（5）：591－598.

[4]Rife D C，Vincent G A.Use of the discrete Fourier transform in the measurement of frequencies and levels of tones[J].Bell syst.Tech.J.，1970，49：197－228.

[5]王旭東，劉渝，鄧振淼.基于修正Rife算法的正弦波頻率估計及FPGA實現[J].系統工程與電子技術，2008，30（4）：621－624.

[6]Altera.FFT MegaCore Function v8.0user Guide[EB／OL].http：／／www.altera.com.

[7]齊國清.離散實正弦信號參數估計的Cramer Rao方差下限[J].數據采集與處理，2003，18（2）：151－155.