唐濟(jì)遠(yuǎn) 袁春姍

（中國船舶集團(tuán)公司七五〇試驗(yàn)場(chǎng)，昆明，650000）

對(duì)于水下探測(cè)設(shè)備來說，信號(hào)檢測(cè)是難點(diǎn)，也是熱門問題。水聲環(huán)境復(fù)雜多變，多數(shù)情況下信噪比條件惡劣，頻率估計(jì)精度差，但實(shí)時(shí)性要求又較高。如何在工程實(shí)踐中精確估計(jì)未知信號(hào)的頻率信息，是工程設(shè)計(jì)人員不斷面臨的挑戰(zhàn)。

在不同的信噪比條件下，不同方法所得到的頻率估計(jì)精度各不相同，但是，無論采用哪種頻率估計(jì)方法，其頻率估計(jì)的均方根誤差（Root Mean Square Error，RMSE）都不會(huì)小于克拉美-羅界[1]（Cramer-Rao bound，CLRB）。

文獻(xiàn)[2]給出了基于Rife算法的頻率估計(jì)方法，在對(duì)輸入信號(hào)進(jìn)行一次FFT運(yùn)算后，利用最大譜線及其左邊或右邊的一根次大譜線進(jìn)行插值來確定真實(shí)頻率位置。該算法只需要一次FFT運(yùn)算，運(yùn)算量小，容易硬件實(shí)現(xiàn)。Rife方法在量化頻率中心區(qū)域的頻率估計(jì)誤差很小，但是在量化頻率點(diǎn)附近的誤差卻較大。而修正Rife算法通過對(duì)輸入信號(hào)進(jìn)行頻譜搬移，使得信號(hào)的頻率始終位于量化頻率的中心區(qū)域，提高了頻率估計(jì)精度。

同時(shí)，F(xiàn)PGA的并行性好，文獻(xiàn)[3]介紹了利用FPGA實(shí)現(xiàn)修正 Rife方法的設(shè)計(jì)思路。本文利用XILINX公司的FFT核[4]，將修正Rife方法在FPGA中進(jìn)行了硬件實(shí)現(xiàn)。FPGA布局布線后的時(shí)序仿真表明，本文設(shè)計(jì)實(shí)現(xiàn)了對(duì)采樣率為200 kHz信號(hào)的實(shí)時(shí)頻率估計(jì)功能。

1 方法介紹

1.1 Rife方法介紹

當(dāng)輸入信號(hào)的頻率不在FFT量化頻率點(diǎn)處時(shí)，直接用 FFT運(yùn)算的最大譜線位置來估計(jì)正弦輸入信號(hào)的頻率，存在量化誤差，其誤差范圍為[-fs/2N，fs/2N]（fs為采樣頻率，N為FFT點(diǎn)數(shù)）。Rife方法利用最大譜線以及與其相鄰的次大譜線進(jìn)行插值估計(jì)信號(hào)頻率。設(shè)在高斯白噪聲環(huán)境下的正弦波信號(hào)為

式中，a、fc、Φ0分別為振幅、頻率和初相，Δt為采樣間隔，N為樣本數(shù)。v(t)是實(shí)部和虛部相互獨(dú)立的、方差為2σ2的零均值高斯白噪聲（實(shí)部、虛部都服從N(0,σ2)分布）。對(duì)x(n)作FFT，最大譜線值記為|X(k0)|，次大譜線值記為。則Rife算法所得到的頻率估計(jì)值[1]為

圖1 Rife方法示意圖

文獻(xiàn)[1]證明當(dāng)輸入噪聲為零時(shí)，Rife算法能夠得到精確的頻率估計(jì)結(jié)果。在適當(dāng)?shù)男旁氡葪l件下，當(dāng)fc位于兩個(gè)離散頻率的中心區(qū)域時(shí)，Rife算法性能較好，其頻率估計(jì)均方根誤差接近CLRB，誤差遠(yuǎn)小于 FFT算法；反之，當(dāng)信噪比較低而且fc接近量化頻率點(diǎn)時(shí)，誤差甚至大于FFT算法。這是因?yàn)楫?dāng)fc位于離散頻率中心區(qū)域時(shí)，|X(k0+r)|與|X(k0)|很接近，這時(shí)采用式（2）的內(nèi)插公式具有較高的精度；但是，當(dāng)fc接近于量化頻率點(diǎn)時(shí)，則|X(k0+r)|與|X(k0)|差別很大，|X(k0+r)|很小，在噪聲的影響下，將產(chǎn)生頻率插值的方向性錯(cuò)誤，因此頻率估計(jì)精度降低。

1.2 修正Rife方法及實(shí)現(xiàn)步驟

將Rife方法的頻率插值記為δ1，有[2]

則有δ1∈[0,1/2]。當(dāng)δ1∈[1/3,1/2]時(shí)，Rife 算法精度很高；當(dāng)δ1∈[0,1/3]時(shí)，估計(jì)能力下降。定義δ1∈[1/3,1/2]為量化頻率中心區(qū)域，當(dāng)由Rife得到的頻率估計(jì)結(jié)果位于此中心區(qū)域時(shí)，將其作為最終輸出；否則，將信號(hào)頻移r/3倍量化頻率。信號(hào)頻移采用與復(fù)指數(shù)相乘，從而實(shí)現(xiàn)信號(hào)單方向的頻移。平移后的信號(hào)為

由于是固定平移r/3倍量化頻率，所以復(fù)指數(shù)的計(jì)算可通過查表來實(shí)現(xiàn)，減少了計(jì)算量。對(duì)平移后的信號(hào)再應(yīng)用 Rife算法進(jìn)行頻率估計(jì)，若所得δ1∈[1/3,1/2]則將此時(shí)估計(jì)出的頻率減去r/3倍量化頻率，即為最終的頻率估計(jì)結(jié)果；若所得δ1∈[0,1/3]，則說明頻率移動(dòng)方向錯(cuò)誤，將式（4）中的r取反，再使用Rife方法進(jìn)行頻率估計(jì)，將所得結(jié)果減去r/3倍量化頻率，即為最終的頻率估計(jì)結(jié)果。其運(yùn)算步驟如下：

（1）對(duì)輸入信號(hào)x(n)做FFT運(yùn)算；

（2）計(jì)算δ1，當(dāng)δ1∈[1/3,1/2]，則用式（2）定義的Rife方法進(jìn)行頻率估計(jì)，結(jié)束運(yùn)算。

（3）當(dāng)δ1∈[0,1/3]，用式（4）對(duì)信號(hào)進(jìn)行平移，對(duì)平移后的信號(hào)做 FFT，重新計(jì)算δ1，若δ1∈[1/3,1/2]，用式（2）估計(jì)頻率，并將其減去r/3倍量化頻率作為最終頻率估計(jì)值，結(jié)束運(yùn)算。

（4）當(dāng)重新計(jì)算的δ1∈[0,1/3]，則將r取反，返回第三步。

通過對(duì)原始信號(hào)的頻移，保證了δ1∈[1/3,1/2]，充分發(fā)揮了Rife算法在量化頻率中心區(qū)域精度高的優(yōu)勢(shì)，使得頻率估計(jì)精度大為改善。

2 FPGA實(shí)現(xiàn)方法

2.1 系統(tǒng)設(shè)計(jì)

根據(jù)1.2節(jié)給出的修正Rife方法，在FPGA中完成硬件實(shí)現(xiàn)。硬件設(shè)計(jì)中為了充分發(fā)揮FPGA并行運(yùn)算的優(yōu)勢(shì)，提高運(yùn)算速度，提高修正Rife方法的實(shí)時(shí)性，將輸入信號(hào)分成并行的3路，第2路左向頻移-1/3量化頻率間隔，第3路右向頻移1/3量化頻率間隔，第1路保持原始信號(hào)不變。每路中采用一個(gè)FFT模塊對(duì)各路信號(hào)進(jìn)行FFT運(yùn)算。采用1.2節(jié)介紹的方法計(jì)算3路通道修正量選取其中最接近1/2的值作為有效頻率估計(jì)通道。如果有效頻率估計(jì)通道為第2路，則最終修正量為；如果有效頻率估計(jì)通道為第3路，則最終修正量為；如果有效頻率估計(jì)通道為第1路，則最終修正量為。根據(jù)最終修正量，結(jié)合式（2），對(duì)信號(hào)真實(shí)頻率進(jìn)行估計(jì)。設(shè)計(jì)中采用XILINX公司Artix系列FPGA器件中的XC7A100T2FTG256芯片[4]，在ISE14.7中進(jìn)行編程、綜合、布局布線和時(shí)序仿真。所得設(shè)計(jì)的RTL圖如圖2所示。

圖2 修正Rife方法FPGA實(shí)現(xiàn)RTL圖

系統(tǒng)設(shè)計(jì)主要由3個(gè)FFT模塊、3個(gè)Rife模塊、2個(gè)ROM模塊和一個(gè)判決模塊組成。3個(gè)FFT模塊用來完成3路并行的FFT運(yùn)算；Rife模塊用來實(shí)現(xiàn)Rife方法以及修正量的計(jì)算；2個(gè)ROM模塊用來存放頻移所需復(fù)指數(shù)的實(shí)部和虛部。對(duì)于采樣頻率固定、FFT分析點(diǎn)數(shù)固定的系統(tǒng)來說，頻移系數(shù)為兩個(gè)固定反向的復(fù)指數(shù)，兩個(gè)ROM模塊由于在固定的信號(hào)通路上，所以與固定指數(shù)相乘即可。判決模塊對(duì) 3路并行計(jì)算的結(jié)果進(jìn)行判決，選取最大值，并相應(yīng)的加減頻率，給出最終的頻率分析結(jié)果。FFT模塊可通過XILINX公司FFT核的參數(shù)來完成設(shè)計(jì)，其他模塊采用VHDL硬件描述語言來實(shí)現(xiàn)。

2.2 硬件設(shè)計(jì)資源消耗分析

在ISE14.7軟件中建立工程，經(jīng)綜合、布局布線后的資源使用情況如表 1所示。采用XC7A100T2FTG256芯片能夠在一片 FPGA內(nèi)實(shí)現(xiàn)并行結(jié)構(gòu)的修正Rife方法。這種并行、流水的處理方式，加快了處理速度，提高了系統(tǒng)的實(shí)時(shí)性，適合在水聲信號(hào)檢測(cè)、水聲探測(cè)、水聲誘餌等信噪比較低、系統(tǒng)實(shí)現(xiàn)性要求高的領(lǐng)域應(yīng)用。該設(shè)計(jì)的缺點(diǎn)是消耗的邏輯單元和DSP模塊比較多。

表1 FPGA主要資源消耗表

2.3 設(shè)計(jì)后仿真

將設(shè)計(jì)在 Modelsim軟件中進(jìn)行布線布局，得到設(shè)計(jì)的延遲模型，再用Matlab軟件生成仿真所需的向量文件（*.vec文件）。系統(tǒng)輸入兩段頻率分別為（65.15/N）×fs和（110.5/N）×fs的信號(hào)，長度都為256，信噪比SNR=-3 dB。data_real_in表示輸入信號(hào)，位寬16 bit。系統(tǒng)時(shí)鐘信號(hào)為clk，復(fù)位信號(hào)為reset。系統(tǒng)輸出信號(hào)為k0、r和deta，其中k0和r的含義與式（2）中的一致，r的高、低電平代表取正、取負(fù)。deta為Rife算法的最終修正量，位寬為16 bit，用來表示16位無符號(hào)定點(diǎn)小數(shù)。當(dāng)修正量為1時(shí)對(duì)應(yīng)deta值為65 535。修正Rife方法的時(shí)序仿真波形如圖3所示，局部信息放大如圖4所示。根據(jù)圖中顯示，第一段信號(hào)k0=65、r=+1、date=8782，則實(shí)際估計(jì)的頻率k0+r×deta/65 535=65.134, 單位是1/N×fs；第二段信號(hào)k0=111、r=-1、date=30 867，則實(shí)際估計(jì)的頻率為k0+r×deta/65 535=110.529, 單位是 1/N×fs。時(shí)序仿真表明，設(shè)計(jì)最高運(yùn)行速度可滿足采樣率200 kHz信號(hào)的處理要求，數(shù)據(jù)無延遲，能夠?qū)崿F(xiàn)并行、流水、實(shí)時(shí)處理，可以較為準(zhǔn)確的估計(jì)真實(shí)信號(hào)頻率。

圖3 修正Rife方法FPGA仿真圖

圖4 修正Rife方法FPGA仿真細(xì)節(jié)圖

3 小結(jié)

本文分析了原始Rife方法和修正Rife方法各自的優(yōu)缺點(diǎn)，說明了修正Rife方法在量化頻率點(diǎn)附近對(duì)頻率估計(jì)精度的提高效果。對(duì)修正Rife方法在FPGA中的硬件設(shè)計(jì)方法進(jìn)行詳細(xì)介紹。充分利用了FPGA器件的優(yōu)勢(shì)結(jié)構(gòu)，對(duì)修正Rife方法進(jìn)行并行運(yùn)算設(shè)計(jì)。最后，給出了FPGA資源消耗情況以及后仿真結(jié)果。該設(shè)計(jì)方法可以應(yīng)用于多種水下信號(hào)檢測(cè)設(shè)備中，在低信噪比的水聲環(huán)境、實(shí)時(shí)性要求高的水聲檢測(cè)系統(tǒng)中有較好頻率估計(jì)效果。