蘇俊亞,張霞
(中國(guó)傳媒大學(xué) 理工學(xué)部,北京 100024)
基于全相位FFT的頻率估計(jì)及其FPGA實(shí)現(xiàn)
蘇俊亞,張霞
(中國(guó)傳媒大學(xué) 理工學(xué)部,北京 100024)
信號(hào)截?cái)鄷?huì)引起頻譜的泄露,導(dǎo)致頻率測(cè)量出現(xiàn)偏差。相對(duì)于傳統(tǒng)FFT,全相位FFT能夠抑制頻譜泄露。本文針對(duì)一種諧振式紅外傳感器,首先對(duì)檢測(cè)信號(hào)進(jìn)行全相位FFT,然后利用傳統(tǒng)FFT的頻譜對(duì)全相位FFT的頻譜進(jìn)行校正,實(shí)現(xiàn)對(duì)頻率的精確估計(jì)。將此測(cè)頻方法在FPGA上進(jìn)行硬件綜合,利用Verilog HDL 語(yǔ)言,設(shè)計(jì)了頻率測(cè)量電路。然后對(duì)所設(shè)計(jì)的頻率測(cè)量電路進(jìn)行仿真,仿真結(jié)果顯示該方法在使用單精度浮點(diǎn)數(shù)進(jìn)行運(yùn)算時(shí)在MHz量級(jí)上可分辨1Hz的信號(hào)頻率。
全相位FFT;相位差;FPGA;Modelsim仿真
頻率測(cè)量是眾多電子技術(shù)中的一種基本技術(shù)[1],在通信、導(dǎo)航系統(tǒng)、空間科學(xué)計(jì)算等領(lǐng)域有著廣泛的應(yīng)用。對(duì)于頻率的測(cè)量,可以將時(shí)域信號(hào)進(jìn)行FFT,然后在頻域?qū)π盘?hào)頻率進(jìn)行估計(jì),原始的信號(hào)序列較長(zhǎng),信號(hào)處理設(shè)備處理的長(zhǎng)度有限,在進(jìn)行FFT變換之前需對(duì)信號(hào)進(jìn)行截?cái)啵瑢?dǎo)致了信號(hào)的頻譜泄露。為了抑制信號(hào)頻譜泄露,通常情況下可以對(duì)信號(hào)進(jìn)行加窗,如漢寧窗、凱澤窗、布萊克曼窗等,這些窗函數(shù)的傅里葉頻譜的第一旁瓣衰減大于矩形窗[2],但是主瓣寬度也大于矩形窗,加窗處理在抑制旁瓣泄露的同時(shí)會(huì)出現(xiàn)主瓣模糊的問(wèn)題,因此加窗處理不能較好地解決頻譜泄露的問(wèn)題,導(dǎo)致頻率測(cè)量的精度不高。
針對(duì)以上問(wèn)題,文獻(xiàn)[3]提出了一種數(shù)字信號(hào)處理的新方法,即全相位快速傅里葉變換(all phase Fast Fourier Transformation,apFFT),這種信號(hào)處理的方法在抑制頻譜泄露的同時(shí),不會(huì)導(dǎo)致主瓣模糊的問(wèn)題。文獻(xiàn)[4]提出了利用傳統(tǒng)FFT的頻譜對(duì)apFFT的頻譜進(jìn)行校正的方法,即FFT/apFFT綜合相位差校正法,該校正法提高了頻率測(cè)量精度。本文針對(duì)一種諧振式紅外傳感器[5],將FFT/apFFT綜合相位差校正法在FPGA上進(jìn)行硬件綜合,采用Verilog HDL 語(yǔ)言,利用開(kāi)發(fā)軟件提供的IP核,設(shè)計(jì)了頻率測(cè)量電路。然后對(duì)所設(shè)計(jì)的頻率測(cè)量電路使用Modelsim進(jìn)行仿真,實(shí)現(xiàn)了對(duì)諧振式紅外傳感器頻率的估計(jì)。
FFT/apFFT綜合相位差校正法[3]是對(duì)信號(hào)分別進(jìn)行全相位FFT和傳統(tǒng)FFT變換。全相位FFT是將采樣的2N-1點(diǎn)數(shù)據(jù)經(jīng)過(guò)預(yù)處理后得到N個(gè)數(shù)據(jù),然后再進(jìn)行傳統(tǒng)的FFT。由于全相位FFT和傳統(tǒng)FFT在齊次、疊加、時(shí)不變、頻移等性質(zhì)存在相似性,因此,可以采用傳統(tǒng)FFT對(duì)全相位FFT的頻譜進(jìn)行校正。
在實(shí)際應(yīng)用中,信號(hào)常為正弦或余弦信號(hào),根據(jù)歐拉公式可知正弦和余弦信號(hào)由復(fù)指數(shù)信號(hào)經(jīng)過(guò)線性疊加而得,因此只需研究復(fù)指數(shù)信號(hào)。一個(gè)幅值為A,初相位為θ的單頻復(fù)指數(shù)信號(hào)可表示為
x(n)=Aej[w*(n-n0)+θ]n∈[0,N-1]
(1)
對(duì)其進(jìn)行傅里葉變換,并對(duì)變換結(jié)果在wk=kΔw=2kπ/N上進(jìn)行采樣,可得其FFT頻譜為
X(k)=Aej(θ-w*n0)2πδ(kΔw-w*)
(2)
窗函數(shù)序列為w(n),n∈[0,N-1],通常窗函數(shù)具有對(duì)稱性,滿足w(n)=w(N-1-n),其FFT譜為
W(k)=Wg(kΔw)e-jτkΔw
(3)
其中τ=(N-1)/2,對(duì)信號(hào)加窗后的FFT譜為
=Aej(θ-w*n0)δ(kΔw-w*)*Wg(kΔw)e-jτkΔw
=Aej(θ-w*n0)Wg[kΔw-w*]e-jτ(kΔw-w*),k∈[0,N-1]
(4)
對(duì)(1)式中的x(n)進(jìn)行前、后窗函數(shù)序列均為w(n)的雙窗全相位FFT,變換后的頻譜為[5]
(5)
假設(shè)最大譜線所對(duì)應(yīng)k值為k*,經(jīng)過(guò)傳統(tǒng)FFT的(4)式與經(jīng)過(guò)全相位FFT的(5)式的相位譜差值為
Δφ=φx(k*)-φy(k*)
=[θ-w*n0-τ(k*Δw-w*)]-(θ-w*n0)
=-τ(k*Δw-w*)
(6)
由上式可得
w*=Δφ/τ+k*Δw,τ=(N-1)/2
(7)
從(7)式可看出只要求出兩種譜線在主譜線上的相位差值Δφ,就可以得出信號(hào)頻率的估計(jì)值。通過(guò)上述對(duì)FFT/apFFT綜合相位差校正法的介紹,可得頻率估計(jì)的流程如圖1所示[6]。
圖1 基于FFT/apFFT綜合相位差校正法的頻率估計(jì)流程圖
此頻率估計(jì)方法包括傳統(tǒng)加窗FFT和雙窗全相位FFT,分別求出兩者主譜線的相位值,其差值除以(N-1)/2可得頻偏值,頻偏值與主譜線的頻率值相加后可得輸入信號(hào)頻率的估計(jì)值。FFT/apFFT綜合相位差法與其他頻譜校正方法相比,其綜合了FFT和apFFT兩者的優(yōu)點(diǎn),不但具有傳統(tǒng)相位差法的所有性能,而且兼?zhèn)淞薬pFFT的頻譜泄露小,可檢測(cè)微弱信號(hào)的優(yōu)點(diǎn),其效率較高,性能較好[7]。
選用Cyclone IV E系列的FPGA芯片不僅具有低成本、低功耗的優(yōu)點(diǎn),而且具有豐富的邏輯單元、I/O、存儲(chǔ)器等資源,足以實(shí)現(xiàn)基于FFT/apFFT綜合相位差校正法的頻率測(cè)量方案;設(shè)計(jì)所用開(kāi)發(fā)軟件為QuartusⅡ14.1,該版本的軟件設(shè)計(jì)了較多兼容Cyclone IV E 系列的IP核,利用這些已有的IP核進(jìn)行設(shè)計(jì)縮短了開(kāi)發(fā)時(shí)間;對(duì)于電路的仿真采用Modelsim SE 10.4,利用該版本可以更加快速實(shí)現(xiàn)系統(tǒng)的功能仿真;數(shù)據(jù)運(yùn)算全部采用IEEE-754標(biāo)準(zhǔn)規(guī)定的單精度浮點(diǎn)數(shù),與定點(diǎn)數(shù)運(yùn)算相比,浮點(diǎn)數(shù)可以表示的數(shù)據(jù)范圍更寬,計(jì)算的精度更高。FPGA開(kāi)發(fā)過(guò)程中依據(jù)自頂而下、模塊化的設(shè)計(jì)理念,圖2是整體設(shè)計(jì)框圖,進(jìn)一步將分析每一模塊電路,并使用Modelsim對(duì)其進(jìn)行仿真,驗(yàn)證是否實(shí)現(xiàn)電路功能。
圖2 FPGA設(shè)計(jì)整體框圖
在全相位FFT數(shù)據(jù)預(yù)處理中,根據(jù)不同卷積窗函數(shù)可分為無(wú)窗全相位數(shù)據(jù)預(yù)處理、單窗全相位數(shù)據(jù)預(yù)處理和雙窗全相位數(shù)據(jù)預(yù)處理,由圖1可知此校正算法采用雙窗全相位數(shù)據(jù)預(yù)處理。從主瓣寬度和旁瓣衰減幅度兩方面考慮,雙窗全相位數(shù)據(jù)預(yù)處理中的前窗函數(shù)和后窗函數(shù)均采用漢寧窗。在數(shù)據(jù)預(yù)處理中,使用MATLAB生成歸一化的卷積窗序列,作為rom的初始化文件,加載到單端口rom中。整個(gè)數(shù)據(jù)預(yù)處理模塊包括rom模塊、乘法模塊、fifo模塊、reg_x0模塊、加法模塊、控制模塊以及數(shù)據(jù)選擇模塊。其中控制模塊為rom模塊提供地址信號(hào),為fifo模塊和reg_x0模塊提供讀寫(xiě)使能信號(hào),為數(shù)據(jù)選擇模塊提供選擇使能信號(hào),是數(shù)據(jù)預(yù)處理電路的核心模塊,協(xié)調(diào)整個(gè)電路工作。圖3是數(shù)據(jù)預(yù)處理模塊的原理圖。
傳統(tǒng)FFT數(shù)據(jù)處理只需要將采樣后的N點(diǎn)數(shù)據(jù)與N點(diǎn)窗函數(shù)序列相乘,相乘后的數(shù)據(jù)進(jìn)行FFT變換即可。與全相位數(shù)據(jù)預(yù)處理相同,窗函數(shù)也選擇漢寧窗,在MATLAB中生成漢寧窗序列,并經(jīng)過(guò)編碼變?yōu)閱尉雀↑c(diǎn)數(shù)的數(shù)據(jù)形式,將其存儲(chǔ)為MIF格式的文件,作為FPGA中rom的初始化文件,加載到單端口rom中,與輸入未知信號(hào)進(jìn)行運(yùn)算。傳統(tǒng)FFT數(shù)據(jù)處理模塊包含三個(gè)模塊,control模塊、rom模塊和mult模塊。control模塊控制rom模塊的地址,地址每變換一次,rom會(huì)輸出一個(gè)窗函數(shù)數(shù)據(jù),rom的輸出會(huì)延遲一個(gè)時(shí)鐘,在testbench中控制未知信號(hào)的輸入時(shí)間,與rom的輸出保持一致,同步送入mult模塊中,實(shí)現(xiàn)對(duì)未知信號(hào)的“加窗”處理。圖4是傳統(tǒng)FFT數(shù)據(jù)處理的原理圖。
FFT變換模塊采用Altera公司研發(fā)的FFT IP核,F(xiàn)FT v14.1可實(shí)現(xiàn)定點(diǎn)、單浮點(diǎn)、塊浮點(diǎn)的FFT和IFFT運(yùn)算[8],變換的長(zhǎng)度為2m(3≤m≤18),數(shù)據(jù)流有流、變量流、緩存突發(fā)、突發(fā)四種結(jié)構(gòu)。本次設(shè)計(jì)采用IEEE-754標(biāo)準(zhǔn)規(guī)定的單精度浮點(diǎn)數(shù),此格式的數(shù)據(jù)位數(shù)為32位,F(xiàn)FT變換的長(zhǎng)度選擇1024點(diǎn),數(shù)據(jù)流結(jié)構(gòu)選擇變量流。在對(duì)FFT核中各參數(shù)分析的基礎(chǔ)上,將其進(jìn)行配置,保存配置后即可自動(dòng)生成Verilog VDL文件,將生成的文件加載到創(chuàng)建的工程中,即可在電路設(shè)計(jì)的過(guò)程中使用元件例化語(yǔ)句對(duì)生成的模塊進(jìn)行調(diào)用實(shí)現(xiàn)其功能。具體配置FFT核的參數(shù)如圖5所示。
求幅值和相位模塊的目的是求出FFT變換后輸出復(fù)數(shù)的幅值和相位,運(yùn)算的過(guò)程均采用單精度浮點(diǎn)數(shù)。在此模塊中調(diào)用Quartus Ⅱ14.1中參數(shù)化模塊庫(kù)(LPM,Library of Parameterized Modules)里面的乘法、加法、除法和反正切模塊,這些模塊對(duì)Altera器件進(jìn)行了優(yōu)化,不耗用器件的邏輯單元,加快了設(shè)計(jì)的開(kāi)發(fā)。求復(fù)數(shù)幅值采用乘法和加法模塊實(shí)現(xiàn)運(yùn)算,求復(fù)數(shù)的相位采用除法和反正切模塊實(shí)現(xiàn)運(yùn)算,反正切模塊中角度的范圍為-π/2~π/2,但是相位差校正算法中使用的角度范圍為-π~π,因此必須采用校正方法對(duì)反正切模塊輸出結(jié)果進(jìn)行校正,在代碼中增加判斷依據(jù),具體的描述如下。
圖3 全相位FFT數(shù)據(jù)預(yù)處理原理圖
圖4 傳統(tǒng)FFT數(shù)據(jù)處理的原理圖
圖5 FFT核的配置
(8)
在求幅值和相位模塊中實(shí)現(xiàn)編號(hào)、幅值、相位三者同步后才能將結(jié)果輸入到下一模塊中進(jìn)行運(yùn)算,設(shè)計(jì)中采用觸發(fā)器作為結(jié)果延遲單元,保證了結(jié)果同步輸出。使用Modelsim對(duì)此模塊進(jìn)行仿真,仿真結(jié)果如圖6所示,從圖中可看出result_num、result_mod、result_phase同步,達(dá)到設(shè)計(jì)要求。
圖6 求幅值和相位模塊的仿真結(jié)果
求最大譜線模塊中,只對(duì)輸入的幅值信號(hào)進(jìn)行比較,幅值信號(hào)是大于等于零的,因此可以設(shè)置一個(gè)初始值為零的臨時(shí)變量與輸入的第一個(gè)數(shù)據(jù)進(jìn)行比較,每一個(gè)CLK輸入一個(gè)數(shù)據(jù),經(jīng)過(guò)比較后將較大的數(shù)據(jù)放入臨時(shí)變量中用于下一個(gè)時(shí)鐘的比較。通過(guò)對(duì)全部幅值進(jìn)行比較,求出最大幅值,并保存最大幅值所對(duì)應(yīng)的相位值和編號(hào)值。求最大幅值的流程圖見(jiàn)圖7。使用Modelsim對(duì)此模塊進(jìn)行仿真,仿真結(jié)果如圖8所示,從圖中可以看出,經(jīng)過(guò)幅值的不斷比較后得到主譜線的編號(hào)值、幅度值和相位值。
圖7 求最大幅值的流程圖
上述各個(gè)模塊連接在一起可組成全相位FFT和傳統(tǒng)FFT求主譜線相位值的電路設(shè)計(jì)。圖9是全相位FFT求主譜線相位值電路設(shè)計(jì)的頂層原理
圖8 求最大譜線模塊的仿真結(jié)果
圖,將未知信號(hào)的2N-1點(diǎn)采樣數(shù)據(jù)送入到如圖所示的電路中,經(jīng)過(guò)數(shù)據(jù)的預(yù)處理,N點(diǎn)FFT變換,求出幅值和相位,然后求出主譜線所對(duì)應(yīng)的幅度值、相位值和編號(hào)值。圖10是傳統(tǒng)FFT求主譜線相位值電路設(shè)計(jì)的頂層原理圖,將未知信號(hào)中包含全相位FFT所需數(shù)據(jù)的第N點(diǎn)到第2N-1點(diǎn)采樣數(shù)據(jù)送入圖中所示電路中,經(jīng)過(guò)“加窗”處理,N點(diǎn)FFT變換,求幅值和相位,然后求出經(jīng)過(guò)普通FFT運(yùn)算后主譜線的幅度值、相位值和編號(hào)值。
在計(jì)算頻率模塊中,全相位FFT和傳統(tǒng)FFT的主譜線相位值和編號(hào)值作為此模塊的輸入數(shù)據(jù),在此模塊中首先計(jì)算傳統(tǒng)FFT和全相位FFT主譜線的相位差值,然后除以(N-1)/2,最后與序列k值和fs/N的乘積進(jìn)行加法運(yùn)算,運(yùn)算的結(jié)果是對(duì)未知信號(hào)頻率的估計(jì)值。假設(shè)傳統(tǒng)FFT主譜線對(duì)應(yīng)的相位角度值為p1,全相位FFT主譜線對(duì)應(yīng)相位角度值為p2,主譜線對(duì)應(yīng)的編號(hào)為k,采樣頻率為fs,變換點(diǎn)數(shù)為N,則對(duì)未知信號(hào)頻率的估計(jì)值為式9[9]。
圖9 全相位FFT求主譜線相位值電路設(shè)計(jì)的頂層原理圖
圖10 傳統(tǒng)FFT求主譜線相位值電路設(shè)計(jì)的頂層原理圖
圖11 計(jì)算頻率模塊的RTL視圖
(9)
計(jì)算頻率模塊的RTL視圖見(jiàn)圖11。圖12是當(dāng)輸入信號(hào)為4.5MHz時(shí),使用Modelsim對(duì)此模塊仿真的結(jié)果。從圖中可以看出,設(shè)計(jì)的電路輸出結(jié)果為未知信號(hào)的頻率,實(shí)現(xiàn)了頻率測(cè)量電路的功能。
圖12 計(jì)算頻率模塊的仿真結(jié)果
本文所測(cè)量的諧振式紅外傳感器的固有頻率約為4.5MHz,因此頻率測(cè)試在4.5MHz左右展開(kāi)。測(cè)試使用信號(hào)為正弦信號(hào),初相位為零,采樣頻率為45MHz,進(jìn)行1024點(diǎn)FFT變換;測(cè)試信號(hào)由MATLAB產(chǎn)生,并在MATLAB中轉(zhuǎn)化成IEEE-754標(biāo)準(zhǔn)的單精度浮點(diǎn)數(shù),寫(xiě)入文本文件作為電路設(shè)計(jì)輸入。使用Modelsim對(duì)不同頻率的輸入信號(hào)進(jìn)行仿真,仿真的結(jié)果如表1所示。
表1 頻率測(cè)量電路仿真結(jié)果
從表中數(shù)據(jù)看出,仿真結(jié)果與實(shí)際信號(hào)頻率相一致。本論文的頻率測(cè)量是針對(duì)參考資料[5]的諧振式紅外傳感器固有頻率進(jìn)行檢測(cè),能夠在MHz量級(jí)分辨1Hz即可滿足要求。顯然,仿真結(jié)果達(dá)到此要求。
本文在FPGA上綜合了FFT/apFFT綜合相位差校正法,在綜合算法過(guò)程中首先進(jìn)行了全相位FFT變換,然后使用傳統(tǒng)FFT對(duì)全相位FFT的頻譜進(jìn)行校正,校正后的頻譜可得到頻率值的精確估計(jì)。在FPGA硬件綜合過(guò)程中使用了Quartus Ⅱ14.1提供的IP核,優(yōu)化了設(shè)計(jì),縮短了開(kāi)發(fā)時(shí)間。FPGA開(kāi)發(fā)是基于自頂而下、模塊化的理念,使用Modelsim軟件對(duì)設(shè)計(jì)的模塊電路進(jìn)行自下而上、層次化遞加的功能仿真,最后對(duì)整個(gè)電路進(jìn)行功能仿真,對(duì)仿真的結(jié)果分析可知,這種校正方法可以實(shí)現(xiàn)對(duì)諧振式紅外傳感器固有頻率的測(cè)量。
[1]陳麗鋒,孟瑞,馮希.頻率測(cè)量研究綜述[J].電力電子,2011(2):50-53.
[2]周新星,王典洪.窗函數(shù)用于頻譜分析的研究[J].軟件時(shí)空,2009,25(12-3):185.
[3]王兆華,黃翔東.數(shù)字信號(hào)全相位譜分析與濾波技術(shù)[M].北京:電子工業(yè)出版社,2009,63-69.
[4]黃翔東,王兆華.基于全相位頻譜分析的相位差校正法[J].電子與信息學(xué)報(bào),2008,30(2):293-297.
[5]Xia Zhang,Dacheng Zhang.Thermal expansion mismatch resonant infrared sensor and array [J]. Optical Engineering,2015,54(8):087102.
[6]王兆華,侯正信,蘇飛.全相位FFT頻譜分析[J].通信學(xué)報(bào),2003,24(11A):6-19.
[7]黃翔東.全相位數(shù)字信號(hào)處理[D]. 天津:天津大學(xué),2006.
[8]劉東華.Altera系列FPGA芯片IP核詳解[M].北京:電子工業(yè)出版社,2014.
[9]任麗棉.基于全相位FFT的高精度頻率計(jì)系統(tǒng)研究[D].天津:天津大學(xué),2009.
(責(zé)任編輯:宋金寶)
FrequencyEstimationBasedonFullPhaseFFTandItsImplementationonFPGA
SU Jun-ya,ZHANG Xia
(Faculty of Science and Technology,Communication University of China,Beijing 100024,China)
The signal cut off will cause a leakage of spectrum,resulting in frequency measurement deviation. Compared to the traditional FFT,all-phase FFT can inhibit the leakage of the spectrum. For a resonant infrared sensor in this paper,firstly,the full-phase FFT is performed for detection signal,and then conventional FFT can be used to correct the spectrum of the all-phase FFT to achieve the accurate estimation of frequency. This frequency measurement method is synthesized on the FPGA. Using the Verilog HDL language,a frequency measurement circuit is designed. Then,the frequency measurement circuit is simulated. By analyzing the simulation results,it is found that the method can distinguish the signal frequency of 1Hz around 4.5 MHz when single-precision floating-point is used to operation.
all-phase FFT;phase difference;FPGA;Modelsim simulation
TM935
A
1673-4793(2017)05-0050-07
2017-06-28
蘇俊亞(1988-),女(漢族),河北邯鄲人,中國(guó)傳媒大學(xué)碩士研究生.E-mail:sujunya@cuc.edu.cn
中國(guó)傳媒大學(xué)學(xué)報(bào)(自然科學(xué)版)2017年5期